From 160d13927a97dc06831157b32f5614c54389cc47 Mon Sep 17 00:00:00 2001 From: ocayon Date: Tue, 3 Nov 2026 13:47:20 +0100 Subject: [PATCH 1/3] perf(utils): allocation-free numba filament kernels (~70x faster AIC) The induced-velocity kernels carried vectors as length-3 arrays; in numba every temporary (np.cross, r1 - r0, np.zeros(3)) is a heap allocation, and those allocations were ~95% of the cost: ~1.1 us per filament evaluation. They now carry 3-tuples, which numba keeps on the stack (~15-19 ns). Formulas, branch structure and summation order are unchanged, so the matrices still match the Filament.py reference to round-off (the equivalence tests at rtol 1e-12 pass). error_model="numpy" keeps NumPy's inf/nan semantics for zero denominators instead of numba's scalar ZeroDivisionError, so degenerate geometry behaves as before. V3 kite, 36 panels: AIC 7.3 -> 0.10 ms; casadi_newton solve 44 -> 24 ms. Affects assemble_AIC_matrices, assemble_curved_bound_correction_AIC and induced_velocity_at_points. Numerical, not physical; outputs unchanged. Generated with Claude Code --- src/VSM/core/utils.py | 363 ++++++++++++++++++++++++------------------ 1 file changed, 207 insertions(+), 156 deletions(-) diff --git a/src/VSM/core/utils.py b/src/VSM/core/utils.py index d4c4643..10c3630 100644 --- a/src/VSM/core/utils.py +++ b/src/VSM/core/utils.py @@ -30,156 +30,210 @@ def jit_dot(a, b): # round-off; Filament.py remains the readable reference implementation. # The only intentional difference: the bound-vortex in-core branch cannot # emit the logging.info diagnostic from compiled code. +# +# Vectors are carried as 3-tuples of floats, not length-3 arrays: numba +# keeps tuples on the stack, whereas every small-array temporary (np.cross, +# r1 - r0, np.zeros(3)) is a heap allocation, and those allocations were +# ~95% of the kernel cost (~1.1 us -> ~25 ns per filament evaluation). +# error_model="numpy" keeps NumPy's inf/nan semantics for a zero +# denominator (numba's default would raise ZeroDivisionError on scalars), +# so degenerate geometry behaves as it did with the array kernels. # ===================================================================== _ALPHA0_OSEEN = 1.25643 # Oseen parameter (Filament._alpha0) _NU_AIR = 1.48e-5 # kinematic viscosity of air [m^2/s] (Filament._nu) +_ZERO3 = (0.0, 0.0, 0.0) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _row(a, i): + """Row ``i`` of an (n, 3) array as a 3-tuple.""" + return (a[i, 0], a[i, 1], a[i, 2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _add(a, b): + return (a[0] + b[0], a[1] + b[1], a[2] + b[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _sub(a, b): + return (a[0] - b[0], a[1] - b[1], a[2] - b[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _scale(s, a): + return (s * a[0], s * a[1], s * a[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _div(a, s): + return (a[0] / s, a[1] / s, a[2] / s) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _dot(a, b): + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] -@jit(nopython=True, cache=True) + +@jit(nopython=True, cache=True, error_model="numpy") +def _cross(a, b): + return ( + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _norm(a): + return np.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") def _vel_bound_vortex(XV1, XV2, XVP, gamma, core_radius_fraction): """Bound-filament induced velocity (Vatistas core), as in Filament.velocity_3D_bound_vortex.""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - r2 = XVP - XV2 + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + r2 = _sub(XVP, XV2) - r1Xr0 = np.cross(r1, r0) - epsilon = core_radius_fraction * np.linalg.norm(r0) - dist = np.linalg.norm(r1Xr0) / np.linalg.norm(r0) + r1Xr0 = _cross(r1, r0) + epsilon = core_radius_fraction * _norm(r0) + dist = _norm(r1Xr0) / _norm(r0) if dist > epsilon: - r1Xr2 = np.cross(r1, r2) - return ( + r1Xr2 = _cross(r1, r2) + return _scale( gamma / (4 * np.pi) - * r1Xr2 - / (np.linalg.norm(r1Xr2) ** 2) - * np.dot(r0, r1 / np.linalg.norm(r1) - r2 / np.linalg.norm(r2)) + / (_norm(r1Xr2) ** 2) + * _dot(r0, _sub(_div(r1, _norm(r1)), _div(r2, _norm(r2)))), + r1Xr2, ) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - r2_radial = r2 - np.dot(r2, r0) * r0 / (np.linalg.norm(r0) ** 2) - r1_proj = np.dot(r1, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r1_radial / np.linalg.norm(r1_radial) - r2_proj = np.dot(r2, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r2_radial / np.linalg.norm(r2_radial) - r1Xr2_proj = np.cross(r1_proj, r2_proj) - vel_ind_proj = ( + r0_sq = _norm(r0) ** 2 + r1_axial = _div(_scale(_dot(r1, r0), r0), r0_sq) + r2_axial = _div(_scale(_dot(r2, r0), r0), r0_sq) + r1_radial = _sub(r1, r1_axial) + r2_radial = _sub(r2, r2_axial) + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r2_proj = _add(r2_axial, _div(_scale(epsilon, r2_radial), _norm(r2_radial))) + r1Xr2_proj = _cross(r1_proj, r2_proj) + vel_ind_proj = _scale( gamma / (4 * np.pi) - * r1Xr2_proj - / (np.linalg.norm(r1Xr2_proj) ** 2) - * np.dot( + / (_norm(r1Xr2_proj) ** 2) + * _dot( r0, - r1_proj / np.linalg.norm(r1_proj) - r2_proj / np.linalg.norm(r2_proj), - ) + _sub(_div(r1_proj, _norm(r1_proj)), _div(r2_proj, _norm(r2_proj))), + ), + r1Xr2_proj, ) - return dist / epsilon * vel_ind_proj + return _scale(dist / epsilon, vel_ind_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_trailing_vortex(XV1, XV2, XVP, gamma, Uinf): """Finite trailing-leg induced velocity (viscous core), as in Filament.velocity_3D_trailing_vortex.""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - r2 = XVP - XV2 + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + r2 = _sub(XVP, XV2) - r_perp = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * np.linalg.norm(r_perp) / Uinf) + r0_sq = _norm(r0) ** 2 + r_perp = _sub(r1, _div(_scale(_dot(r1, r0), r0), r0_sq)) + epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * _norm(r_perp) / Uinf) - r1Xr0 = np.cross(r1, r0) - dist = np.linalg.norm(r1Xr0) / np.linalg.norm(r0) + r1Xr0 = _cross(r1, r0) + dist = _norm(r1Xr0) / _norm(r0) if dist > epsilon: - r1Xr2 = np.cross(r1, r2) - return ( + r1Xr2 = _cross(r1, r2) + return _scale( gamma / (4 * np.pi) - * r1Xr2 - / (np.linalg.norm(r1Xr2) ** 2) - * np.dot(r0, r1 / np.linalg.norm(r1) - r2 / np.linalg.norm(r2)) + / (_norm(r1Xr2) ** 2) + * _dot(r0, _sub(_div(r1, _norm(r1)), _div(r2, _norm(r2)))), + r1Xr2, ) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - r2_radial = r2 - np.dot(r2, r0) * r0 / (np.linalg.norm(r0) ** 2) - r1_proj = np.dot(r1, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r1_radial / np.linalg.norm(r1_radial) - r2_proj = np.dot(r2, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r2_radial / np.linalg.norm(r2_radial) - r1Xr2_proj = np.cross(r1_proj, r2_proj) - vel_ind_proj = ( + r1_axial = _div(_scale(_dot(r1, r0), r0), r0_sq) + r2_axial = _div(_scale(_dot(r2, r0), r0), r0_sq) + r1_radial = _sub(r1, r1_axial) + r2_radial = _sub(r2, r2_axial) + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r2_proj = _add(r2_axial, _div(_scale(epsilon, r2_radial), _norm(r2_radial))) + r1Xr2_proj = _cross(r1_proj, r2_proj) + vel_ind_proj = _scale( gamma / (4 * np.pi) - * r1Xr2_proj - / (np.linalg.norm(r1Xr2_proj) ** 2) - * np.dot( + / (_norm(r1Xr2_proj) ** 2) + * _dot( r0, - r1_proj / np.linalg.norm(r1_proj) - r2_proj / np.linalg.norm(r2_proj), - ) + _sub(_div(r1_proj, _norm(r1_proj)), _div(r2_proj, _norm(r2_proj))), + ), + r1Xr2_proj, ) - return dist / epsilon * vel_ind_proj + return _scale(dist / epsilon, vel_ind_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_semiinfinite(XV1, Vf, XVP, gamma, Uinf, filament_direction): """Semi-infinite trailing-vortex induced velocity, as in Filament.velocity_3D_trailing_vortex_semiinfinite.""" GAMMA = -gamma * filament_direction - r1 = XVP - XV1 - r1XVf = np.cross(r1, Vf) - r_perp = r1 - np.dot(r1, Vf) * Vf + r1 = _sub(XVP, XV1) + r1XVf = _cross(r1, Vf) + r1_axial = _scale(_dot(r1, Vf), Vf) + r_perp = _sub(r1, r1_axial) - if np.linalg.norm(r_perp) < 1e-16: - return np.zeros(3) + if _norm(r_perp) < 1e-16: + return _ZERO3 - epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * np.linalg.norm(r_perp) / Uinf) + epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * _norm(r_perp) / Uinf) - dist = np.linalg.norm(r1XVf) / np.linalg.norm(Vf) + dist = _norm(r1XVf) / _norm(Vf) if dist > epsilon: K = ( GAMMA / 4 / np.pi - / np.linalg.norm(r1XVf) ** 2 - * (1 + np.dot(r1, Vf) / np.linalg.norm(r1)) + / _norm(r1XVf) ** 2 + * (1 + _dot(r1, Vf) / _norm(r1)) ) - return K * r1XVf + return _scale(K, r1XVf) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, Vf) * Vf - r1_proj = np.dot(r1, Vf) * Vf + epsilon * r1_radial / np.linalg.norm(r1_radial) - r1XVf_proj = np.cross(r1_proj, Vf) + r1_radial = r_perp + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r1XVf_proj = _cross(r1_proj, Vf) K = ( GAMMA / 4 / np.pi - / np.linalg.norm(r1XVf_proj) ** 2 - * (1 + np.dot(r1_proj, Vf) / np.linalg.norm(r1_proj)) + / _norm(r1XVf_proj) ** 2 + * (1 + _dot(r1_proj, Vf) / _norm(r1_proj)) ) - return K * r1XVf_proj + return _scale(K, r1XVf_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_bound_2D(XVP, bound_point_1, bound_point_2): """2D bound-vortex correction at the evaluation point, as in Panel.compute_velocity_induced_bound_2D (gamma = 1).""" - r3 = XVP - (bound_point_1 + bound_point_2) / 2 - r0 = bound_point_1 - bound_point_2 - r0_unit_X_r3 = np.cross(r0 / np.linalg.norm(r0), r3) - return (1 / (2 * np.pi)) * (r0_unit_X_r3 / (np.linalg.norm(r0_unit_X_r3) ** 2)) + r3 = _sub(XVP, _div(_add(bound_point_1, bound_point_2), 2.0)) + r0 = _sub(bound_point_1, bound_point_2) + r0_unit_X_r3 = _cross(_div(r0, _norm(r0)), r3) + return _scale(1 / (2 * np.pi), _div(r0_unit_X_r3, _norm(r0_unit_X_r3) ** 2)) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def assemble_AIC_matrices( eval_points, bound_point_1, @@ -212,43 +266,43 @@ def assemble_AIC_matrices( n = eval_points.shape[0] AIC = np.empty((3, n, n)) for icp in range(n): - ep = eval_points[icp] + ep = _row(eval_points, icp) for jring in range(n): + bp1 = _row(bound_point_1, jring) + bp2 = _row(bound_point_2, jring) + te1 = _row(TE_point_1, jring) + te2 = _row(TE_point_2, jring) if evaluation_point_on_bound: - velocity_induced = np.zeros(3) + velocity_induced = _ZERO3 else: velocity_induced = _vel_bound_vortex( - bound_point_2[jring], - bound_point_1[jring], - ep, - 1.0, - core_radius_fraction, + bp2, bp1, ep, 1.0, core_radius_fraction ) - wake_unit = wake_units[jring] + wake_unit = _row(wake_units, jring) wake_speed = wake_speeds[jring] - velocity_induced = velocity_induced + _vel_trailing_vortex( - bound_point_1[jring], TE_point_1[jring], ep, 1.0, wake_speed + velocity_induced = _add( + velocity_induced, _vel_trailing_vortex(bp1, te1, ep, 1.0, wake_speed) ) - velocity_induced = velocity_induced + _vel_trailing_vortex( - TE_point_2[jring], bound_point_2[jring], ep, 1.0, wake_speed + velocity_induced = _add( + velocity_induced, _vel_trailing_vortex(te2, bp2, ep, 1.0, wake_speed) ) - velocity_induced = velocity_induced + _vel_semiinfinite( - TE_point_1[jring], wake_unit, ep, 1.0, wake_speed, 1.0 + velocity_induced = _add( + velocity_induced, + _vel_semiinfinite(te1, wake_unit, ep, 1.0, wake_speed, 1.0), ) - velocity_induced = velocity_induced + _vel_semiinfinite( - TE_point_2[jring], wake_unit, ep, 1.0, wake_speed, -1.0 + velocity_induced = _add( + velocity_induced, + _vel_semiinfinite(te2, wake_unit, ep, 1.0, wake_speed, -1.0), ) if subtract_bound_2D and icp == jring: - velocity_induced = velocity_induced - _vel_bound_2D( - ep, bound_point_1[jring], bound_point_2[jring] - ) + velocity_induced = _sub(velocity_induced, _vel_bound_2D(ep, bp1, bp2)) AIC[0, icp, jring] = velocity_induced[0] AIC[1, icp, jring] = velocity_induced[1] AIC[2, icp, jring] = velocity_induced[2] return AIC -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def assemble_curved_bound_correction_AIC( eval_points, bound_point_1, @@ -315,50 +369,53 @@ def assemble_curved_bound_correction_AIC( for j in range(n): s_node[j + 1] = s_node[j] + L[j] for i in range(n): - ep = eval_points[i] - t = seg[i] / L[i] - c_i = chord_dir[i] + ep = _row(eval_points, i) + t = _div(_row(seg, i), L[i]) + c_i = _row(chord_dir, i) s_P = 0.5 * (s_node[i] + s_node[i + 1]) U = va_norm[i] for j in range(n): - ra = ep + sgn * (s_node[j] - s_P) * t - rb = ep + sgn * (s_node[j + 1] - s_P) * t + ra = _add(ep, _scale(sgn * (s_node[j] - s_P), t)) + rb = _add(ep, _scale(sgn * (s_node[j + 1] - s_P), t)) if fwd: rb2 = ra rb1 = rb else: rb2 = rb rb1 = ra - leg = (1.0 - ac_fraction) * chords[j] * c_i - rte1 = rb1 + leg - rte2 = rb2 + leg - bp1 = bound_point_1[j] - bp2 = bound_point_2[j] - te1 = TE_point_1[j] - te2 = TE_point_2[j] + leg = _scale((1.0 - ac_fraction) * chords[j], c_i) + rte1 = _add(rb1, leg) + rte2 = _add(rb2, leg) + bp1 = _row(bound_point_1, j) + bp2 = _row(bound_point_2, j) + te1 = _row(TE_point_1, j) + te2 = _row(TE_point_2, j) # unrolled legs (single filaments) minus actual legs (single filaments) - v = _vel_trailing_vortex(rb1, rte1, ep, 1.0, U) + _vel_trailing_vortex( - rte2, rb2, ep, 1.0, U + v = _add( + _vel_trailing_vortex(rb1, rte1, ep, 1.0, U), + _vel_trailing_vortex(rte2, rb2, ep, 1.0, U), ) - v = v - _vel_trailing_vortex(bp1, te1, ep, 1.0, U) - _vel_trailing_vortex( - te2, bp2, ep, 1.0, U + v = _sub( + _sub(v, _vel_trailing_vortex(bp1, te1, ep, 1.0, U)), + _vel_trailing_vortex(te2, bp2, ep, 1.0, U), ) # chordwise-spread actual system minus chordwise-spread unrolled system for q in range(n_strips): f = (strip_x[q] - ac_fraction) / (1.0 - ac_fraction) w = strip_w[q] - n1 = bp1 + f * (te1 - bp1) - n2 = bp2 + f * (te2 - bp2) - u1 = rb1 + f * (rte1 - rb1) - u2 = rb2 + f * (rte2 - rb2) - v = v + w * ( - _vel_bound_vortex(n2, n1, ep, 1.0, core_radius_fraction) - - _vel_bound_vortex(u2, u1, ep, 1.0, core_radius_fraction) - + _vel_trailing_vortex(n1, te1, ep, 1.0, U) - + _vel_trailing_vortex(te2, n2, ep, 1.0, U) - - _vel_trailing_vortex(u1, rte1, ep, 1.0, U) - - _vel_trailing_vortex(rte2, u2, ep, 1.0, U) + n1 = _add(bp1, _scale(f, _sub(te1, bp1))) + n2 = _add(bp2, _scale(f, _sub(te2, bp2))) + u1 = _add(rb1, _scale(f, _sub(rte1, rb1))) + u2 = _add(rb2, _scale(f, _sub(rte2, rb2))) + strip = _sub( + _vel_bound_vortex(n2, n1, ep, 1.0, core_radius_fraction), + _vel_bound_vortex(u2, u1, ep, 1.0, core_radius_fraction), ) + strip = _add(strip, _vel_trailing_vortex(n1, te1, ep, 1.0, U)) + strip = _add(strip, _vel_trailing_vortex(te2, n2, ep, 1.0, U)) + strip = _sub(strip, _vel_trailing_vortex(u1, rte1, ep, 1.0, U)) + strip = _sub(strip, _vel_trailing_vortex(rte2, u2, ep, 1.0, U)) + v = _add(v, _scale(w, strip)) A[0, i, j] = v[0] A[1, i, j] = v[1] A[2, i, j] = v[2] @@ -378,18 +435,18 @@ def thin_airfoil_strips(n_strips=24): return np.ascontiguousarray(x, dtype=float), np.ascontiguousarray(w / w.sum(), dtype=float) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _on_segment_line(XV1, XV2, XVP): """True if XVP lies on the line through XV1 and XV2 (to round-off).""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - len0 = np.linalg.norm(r0) + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + len0 = _norm(r0) if len0 == 0.0: return True - return np.linalg.norm(np.cross(r1, r0)) <= 1e-10 * len0 * max(np.linalg.norm(r1), len0) + return _norm(_cross(r1, r0)) <= 1e-10 * len0 * max(_norm(r1), len0) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def induced_velocity_at_points( points, bound_point_1, @@ -413,31 +470,25 @@ def induced_velocity_at_points( n = bound_point_1.shape[0] out = np.zeros((m, 3)) for ip in range(m): - ep = points[ip] - vel = np.zeros(3) + ep = _row(points, ip) + vel = _ZERO3 for j in range(n): g = gamma[j] - wake_unit = wake_units[j] + bp1 = _row(bound_point_1, j) + bp2 = _row(bound_point_2, j) + te1 = _row(TE_point_1, j) + te2 = _row(TE_point_2, j) + wake_unit = _row(wake_units, j) wake_speed = wake_speeds[j] - vel = vel + _vel_bound_vortex( - bound_point_2[j], bound_point_1[j], ep, g, core_radius_fraction - ) + vel = _add(vel, _vel_bound_vortex(bp2, bp1, ep, g, core_radius_fraction)) # a straight filament induces nothing on its own line; skip the # legs the point lies on (the kernel would divide by zero there) - if not _on_segment_line(bound_point_1[j], TE_point_1[j], ep): - vel = vel + _vel_trailing_vortex( - bound_point_1[j], TE_point_1[j], ep, g, wake_speed - ) - if not _on_segment_line(TE_point_2[j], bound_point_2[j], ep): - vel = vel + _vel_trailing_vortex( - TE_point_2[j], bound_point_2[j], ep, g, wake_speed - ) - vel = vel + _vel_semiinfinite( - TE_point_1[j], wake_unit, ep, g, wake_speed, 1.0 - ) - vel = vel + _vel_semiinfinite( - TE_point_2[j], wake_unit, ep, g, wake_speed, -1.0 - ) + if not _on_segment_line(bp1, te1, ep): + vel = _add(vel, _vel_trailing_vortex(bp1, te1, ep, g, wake_speed)) + if not _on_segment_line(te2, bp2, ep): + vel = _add(vel, _vel_trailing_vortex(te2, bp2, ep, g, wake_speed)) + vel = _add(vel, _vel_semiinfinite(te1, wake_unit, ep, g, wake_speed, 1.0)) + vel = _add(vel, _vel_semiinfinite(te2, wake_unit, ep, g, wake_speed, -1.0)) out[ip, 0] = vel[0] out[ip, 1] = vel[1] out[ip, 2] = vel[2] From 424d15d813794bf95ca13936c39f1b2ecbc6ce3d Mon Sep 17 00:00:00 2001 From: ocayon Date: Tue, 3 Nov 2026 14:00:52 +0100 Subject: [PATCH 2/3] feat(core): quarter-chord force law as xp kernels with a public API The section force law moves into VSM.core.panel_kernels, pure functions over an xp namespace (numpy or casadi), so CasADi callers (AWETrim's CasADi trim) use VSM's formula instead of restating it. compute_results calls the same kernels for its per-panel force, so the formula lives once. - BodyAerodynamics.compute_quarter_chord_AIC_matrices(core_radius_fraction, va_norm_array, va_unit_array): public quarter-chord induction, AIC_LLT + the cached curved bound-vortex correction. update_effective_angle_of_ attack_if_VSM now uses it (one source). - panel_kernels: relative_flow, quarter_chord_velocity, direction_angle, section_force_directions, section_force_magnitudes, section_force_vectors, panel_forces, panel_moment_local, attached_trailed_forces, row helpers. compute_results and compute_attached_trailed_vortex_forces call them. - Solver.casadi_panel_force_function(): CasADi function (v_rel, v_qc, axes, chord, width, rho) -> per-panel force, local moment, alpha, alpha_dir, umag, cl, cd, cm, on the same polar lookups as the casadi_newton residual (factored into _casadi_polar_lookups; the residual function is unchanged). Refactor plus API: no physics change. Every numeric solve output matches VSM 1201f6f to <= 1.5e-14 relative for is_aoa_corrected True and False on a straight, a swept, an arc wing and the LEI V3 kite (frozen reference in tests/BodyAerodynamics/data, written with 1201f6f). CasADi kernels equal the NumPy ones to 1e-12; a CasADi graph with the induction matrices as parameters reproduces compute_results panel forces to 2e-16. Generated with Claude Code --- src/VSM/core/BodyAerodynamics.py | 194 ++++++----- src/VSM/core/Solver.py | 215 +++++++++---- src/VSM/core/panel_kernels.py | 236 ++++++++++++++ .../data/panel_kernel_reference.npz | Bin 0 -> 215246 bytes tests/BodyAerodynamics/panel_kernel_cases.py | 115 +++++++ tests/BodyAerodynamics/test_panel_kernels.py | 303 ++++++++++++++++++ 6 files changed, 904 insertions(+), 159 deletions(-) create mode 100644 src/VSM/core/panel_kernels.py create mode 100644 tests/BodyAerodynamics/data/panel_kernel_reference.npz create mode 100644 tests/BodyAerodynamics/panel_kernel_cases.py create mode 100644 tests/BodyAerodynamics/test_panel_kernels.py diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index 4eaf8ac..467cc95 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -17,6 +17,7 @@ thin_airfoil_strips, ) from . import jit_cross, jit_norm, jit_dot +from VSM.core import panel_kernels def _rigid_rotation_between(points_0, points_1, dirs_0, dirs_1, rel_tol=1e-9): @@ -1331,11 +1332,7 @@ def compute_attached_trailed_vortex_forces( (forces, points): ``(n+1, 3)`` force vectors and the points they act at. """ panels = self.panels - n = len(panels) gamma = np.asarray(gamma, dtype=float).ravel() - gamma_net = np.zeros(n + 1) - gamma_net[:n] += gamma # leg bound_1 -> TE_1 of panel j sits on boundary j - gamma_net[1:] -= gamma # leg TE_2 -> bound_2 of panel j sits on boundary j+1 bound_point_1 = np.ascontiguousarray( [p.bound_point_1 for p in panels], dtype=float @@ -1355,13 +1352,7 @@ def compute_attached_trailed_vortex_forces( t = (cp - ac) / (1.0 - ac) points = np.ascontiguousarray(bound + t * legs) - # inflow at a boundary: mean of the neighbouring panels' inflow va_array = np.asarray(va_array, dtype=float) - va_boundary = np.zeros((n + 1, 3)) - va_boundary[:n] += va_array - va_boundary[1:] += va_array - va_boundary[1:n] *= 0.5 - va_speeds = np.linalg.norm(va_array, axis=1) wake_units, wake_speeds = self._wake_directions( va_speeds, va_array / va_speeds[:, None] @@ -1377,8 +1368,9 @@ def compute_attached_trailed_vortex_forces( wake_speeds, float(core_radius_fraction), ) - v_rel = va_boundary + v_ind - forces = rho * gamma_net[:, None] * np.cross(v_rel, legs) + # net leg circulation and boundary inflow (mean of the neighbouring + # panels') are formed in the kernel + forces = panel_kernels.attached_trailed_forces(gamma, va_array, v_ind, legs, rho) return forces, points def compute_results( @@ -1416,10 +1408,6 @@ def compute_results( cl_array[icp] = panel_i.compute_cl(alpha_array[icp]) cd_array[icp], cm_array[icp] = panel_i.compute_cd_cm(alpha_array[icp]) panel_width_array[icp] = panel_i.width - lift = (cl_array * 0.5 * rho * Umag_array**2 * chord_array)[:, np.newaxis] - drag = (cd_array * 0.5 * rho * Umag_array**2 * chord_array)[:, np.newaxis] - moment = (cm_array * 0.5 * rho * Umag_array**2 * chord_array**2)[:, np.newaxis] - # Quarter-chord force directions (TAT3) apply to the VSM model, whose # circulation loop evaluates at 3/4 chord. The LLT model already # evaluates on the quarter chord without any bound-vortex influence @@ -1440,6 +1428,28 @@ def compute_results( alpha_corrected = alpha_array[:, np.newaxis] alpha_uncorrected = alpha_array[:, np.newaxis] + # Section force law (VSM.core.panel_kernels): magnitudes from the polar + # at the 3/4-chord angle and section speed, directions from + # alpha_corrected (the quarter-chord angle, or the 3/4-chord one). + lift_2D_array, drag_2D_array, moment_2D_array, dir_flow_array = ( + panel_kernels.section_force_vectors( + alpha_array, + Umag_array, + cl_array, + cd_array, + cm_array, + chord_array, + x_airf_array, + y_airf_array, + z_airf_array, + rho, + alpha_dir=np.ravel(alpha_corrected), + ) + ) + moment_local_3D_array = panel_kernels.panel_moment_local( + moment_2D_array, panel_width_array, z_airf_array + ) + # Initializing variables cl_prescribed_va_list = [] cd_prescribed_va_list = [] @@ -1496,45 +1506,19 @@ def compute_results( ### Defining panel_variables # Defining directions of airfoil that defines current panel_i z_airf_span = panel_i.z_airf # along the span - y_airf_chord = panel_i.y_airf # along the chord - x_airf_normal_to_chord = panel_i.x_airf # normal to the chord - # TODO: implement these - alpha_corrected_i = alpha_corrected[i] panel_chord = panel_i.chord panel_width = panel_i.width panel_area = panel_chord * panel_width area_all_panels += panel_area - ### Calculate the direction of the induced apparent wind speed to the airfoil orientation - # this is done using the CORRECTED CALCULATED (comes from gamma distribution) angle of attack - # For VSM the correction is applied, and it is the angle of attack, from calculating induced velocities at the 1/4c aerodynamic center location - # For LTT the correction is NOT applied, and it is the angle of attack, from calculating induced velocities at the 3/4c control point - induced_va_airfoil = ( - np.cos(alpha_corrected_i) * y_airf_chord - + np.sin(alpha_corrected_i) * x_airf_normal_to_chord - ) - dir_induced_va_airfoil = induced_va_airfoil / jit_norm(induced_va_airfoil) - - ### Calculate the direction of the lift and drag vectors - # lift is perpendical/normal to induced apparent wind speed - dir_lift_induced_va = jit_cross(dir_induced_va_airfoil, z_airf_span) - dir_lift_induced_va = dir_lift_induced_va / jit_norm(dir_lift_induced_va) - # drag is parallel/tangential to induced apparent wind speed - # (TAT3: the local inner flow direction, which lies in the panel's - # span-perpendicular plane; the global span axis is NOT used here) - dir_drag_induced_va = dir_induced_va_airfoil - - ### Calculating the MAGNITUDE of the lift and drag - # The VSM and LTT methods do NOT differ here, both use the uncorrected angle of attack - # i.e. evaluate the magnitude at the (3/4c) control point - # 2D AIRFOIL aerodynamic forces, so multiplied by chord - lift_induced_va_mag = lift[i] - drag_induced_va_mag = drag[i] - - # panel force VECTOR NORMAL to CALCULATED induced velocity - lift_induced_va = lift_induced_va_mag * dir_lift_induced_va - # panel force VECTOR TANGENTIAL to CALCULATED induced velocity - drag_induced_va = drag_induced_va_mag * dir_drag_induced_va + ### Section force vectors (per unit span) from the force law above: + # lift normal and drag along the direction-defining flow + # dir_induced_va_airfoil (TAT3: the local inner flow direction, in + # the panel's span-perpendicular plane; the global span axis is NOT + # used), magnitudes at the 3/4c control point for VSM and LLT alike. + dir_induced_va_airfoil = dir_flow_array[i] + lift_induced_va = lift_2D_array[i] + drag_induced_va = drag_2D_array[i] ftotal_induced_va = lift_induced_va + drag_induced_va ### Converting forces to prescribed wing va @@ -1689,13 +1673,10 @@ def compute_results( # (1) Panel aerodynamic center in global frame: panel_ac_global = panel_i.aerodynamic_center # 3D [x, y, z] - # (2) Convert local (2D) pitching moment to a 3D vector in global coords. - # Use the axis around which the moment is defined, - # which is the z-axis pointing "spanwise" - moment_axis_global = panel_i.z_airf - - # Scale by panel width if your 'moment[i]' is 2D moment-per-unit-span: - M_local_3D = moment[i] * moment_axis_global * panel_i.width + # (2) Local (2D) pitching moment as a 3D vector in global coords: + # along the z-axis pointing "spanwise", times the width + # (panel_kernels.panel_moment_local). + M_local_3D = moment_local_3D_array[i] # (3) Force application point used for r x F: # use aerodynamic center (1/4c) for all models. @@ -1937,7 +1918,7 @@ def compute_results( ( "cm_panel_dist", np.divide( - moment[:, 0], + moment_2D_array, 0.5 * rho * Umag_array**2 * chord_array**2, out=np.zeros_like(chord_array), where=(Umag_array > 0.0) & (chord_array > 0.0), @@ -2144,39 +2125,39 @@ def _curved_bound_correction_cached( } return AIC_corr - def update_effective_angle_of_attack_if_VSM( - self, - gamma, - core_radius_fraction, - x_airf_array, - y_airf_array, - va_array, - va_norm_array, - va_unit_array, + def compute_quarter_chord_AIC_matrices( + self, core_radius_fraction, va_norm_array, va_unit_array ): - """Updates the angle of attack at the aerodynamic center of each panel, - Calculated at the AERODYNAMIC CENTER, which needs an update for VSM - And can just use the old value for the LLT + """Induction matrices of the force-direction flow at the panel + quarter chords (aerodynamic centres): ``v_qc = va + AIC_qc gamma``. + + The direction-defining flow is the undisturbed relative flow at the + quarter chord (TAT3 in Gaunaa, Li & Pirrung, TORQUE 2026). The LLT + AIC holds the freestream-independent trailed system (legs + wake) + evaluated on the quarter chord, which is Prandtl's exact downwash for + a straight line but misses the influence of a curved or swept bound + vortex. Evaluating the other panels' bound segments on the line itself + is logarithmically singular (Li et al. 2020, "LL-P&S"), so the + curvature/sweep/skew influence is added as the chord-regularised + difference between the actual and the locally unrolled vortex + systems; see assemble_curved_bound_correction_AIC. + + The matrices depend on the geometry and the wake (``va_norm_array``, + ``va_unit_array``, as in :meth:`compute_AIC_matrices`), not on the + circulation, so a caller can freeze them numerically, e.g. as + parameters of a CasADi graph (see ``VSM.core.panel_kernels``). Args: - None + core_radius_fraction (float): Vortex core radius fraction. + va_norm_array (np.ndarray): (n,) panel inflow speeds [m/s]. + va_unit_array (np.ndarray): (n, 3) panel inflow directions. Returns: - None + Tuple[np.ndarray, np.ndarray, np.ndarray]: the x, y and z + components, each (n, n): ``AIC_LLT + AIC_curved_correction``. """ - # The direction-defining flow is the undisturbed relative flow at the - # quarter chord (TAT3 in Gaunaa, Li & Pirrung, TORQUE 2026). The LLT - # AIC holds the freestream-independent trailed system (legs + wake) - # evaluated on the quarter chord, which is Prandtl's exact downwash for - # a straight line but misses the influence of a curved or swept bound - # vortex. Evaluating the other panels' bound segments on the line itself - # is logarithmically singular (Li et al. 2020, "LL-P&S"), so the - # curvature/sweep/skew influence is added as the chord-regularised - # difference between the actual and the locally unrolled vortex - # systems; see assemble_curved_bound_correction_AIC. - aerodynamic_model_type = "LLT" AIC_x, AIC_y, AIC_z = self.compute_AIC_matrices( - aerodynamic_model_type, core_radius_fraction, va_norm_array, va_unit_array + "LLT", core_radius_fraction, va_norm_array, va_unit_array ) eval_points = np.ascontiguousarray( [panel.aerodynamic_center for panel in self.panels], dtype=float @@ -2194,7 +2175,9 @@ def update_effective_angle_of_attack_if_VSM( [panel.TE_point_2 for panel in self.panels], dtype=float ) chords = np.ascontiguousarray([panel.chord for panel in self.panels], dtype=float) - chord_dir = np.ascontiguousarray(y_airf_array, dtype=float) + chord_dir = np.ascontiguousarray( + [panel.y_airf for panel in self.panels], dtype=float + ) # The correction depends on the panel geometry only (the wake terms # cancel between the actual and the unrolled systems, and the leg core # radius has no measurable effect), so it is computed once per shape @@ -2209,21 +2192,34 @@ def update_effective_angle_of_attack_if_VSM( chord_dir, chords, np.ascontiguousarray(va_norm_array, dtype=float).ravel(), float(core_radius_fraction), ) - AIC_x = AIC_x + AIC_corr[0] - AIC_y = AIC_y + AIC_corr[1] - AIC_z = AIC_z + AIC_corr[2] - induced_velocity_all = np.array( - [ - np.matmul(AIC_x, gamma), - np.matmul(AIC_y, gamma), - np.matmul(AIC_z, gamma), - ] - ).T - relative_velocity_array = va_array + induced_velocity_all - v_normal_array = np.sum(x_airf_array * relative_velocity_array, axis=1) - v_tangential_array = np.sum(y_airf_array * relative_velocity_array, axis=1) - alpha_array = np.arctan2(v_normal_array, v_tangential_array) + return AIC_x + AIC_corr[0], AIC_y + AIC_corr[1], AIC_z + AIC_corr[2] + def update_effective_angle_of_attack_if_VSM( + self, + gamma, + core_radius_fraction, + x_airf_array, + y_airf_array, + va_array, + va_norm_array, + va_unit_array, + ): + """Force-direction angle of each panel: the angle of the relative flow + at the quarter chord, ``va + AIC_qc gamma`` with the matrices of + :meth:`compute_quarter_chord_AIC_matrices`. + + Returns: + np.ndarray: (n, 1) angles [rad]. + """ + AIC_x, AIC_y, AIC_z = self.compute_quarter_chord_AIC_matrices( + core_radius_fraction, va_norm_array, va_unit_array + ) + relative_velocity_array = panel_kernels.quarter_chord_velocity( + va_array, gamma, AIC_x, AIC_y, AIC_z + ) + alpha_array = panel_kernels.direction_angle( + relative_velocity_array, x_airf_array, y_airf_array + ) return alpha_array[:, np.newaxis] def compute_line_aerodynamic_force( diff --git a/src/VSM/core/Solver.py b/src/VSM/core/Solver.py index 576f3ba..42d4342 100644 --- a/src/VSM/core/Solver.py +++ b/src/VSM/core/Solver.py @@ -3,6 +3,7 @@ import logging from scipy.linalg import solve_banded from . import jit_cross +from VSM.core import panel_kernels class Solver: @@ -1028,56 +1029,12 @@ def _casadi_newton_key(self, alpha_tables, cl_tables, cd_tables, cm_tables) -> t digest.hexdigest(), ) - def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): - """Build the CasADi function of the PER-PANEL section physics that - :meth:`gamma_loop_casadi_newton` assembles its residual and exact - Jacobian from. - - The circulation residual whose root the Newton loop finds is - - R(gamma) = (I - diag(mu(alpha)) L) gamma - G_raw(gamma), - - exactly the fixed point the ``base`` and ``anderson`` loops iterate - towards (``gamma = (I - diag(mu) L)^-1 G_raw(gamma)``, see - :meth:`_regularize_gamma_target`), including the Li/Gaunaa viscosity - ``mu_i = max(0, -k S Cl'_i / dz_i^2)`` gated on any panel being past - its stall onset (``mu = 0`` without artificial viscosity). Everything - that is not the dense induction ``v_ind = AIC gamma`` is local to a - panel: ``G_raw_i``, ``mu_i`` and ``alpha_i`` depend on ``gamma`` only - through the panel's own relative velocity ``v_rel_i``. So the symbolic - function takes ``v_rel`` (and ``L gamma``) as INPUTS and returns - - h_i = G_raw_i + mu_i (L gamma)_i, c_i = dh_i / dv_rel_i (3 values) - - plus ``mu``, ``alpha`` and ``|v_rel x z|``; the Newton loop then forms - ``R = gamma - h`` and, by the chain rule through ``v_rel = va + AIC - gamma`` (with ``dh_i/d(L gamma)_i = mu_i``), - - J = I - diag(mu) L - sum_k diag(c_k) AIC_k - - in NumPy. This keeps the CasADi graph O(n) -- independent of the - n x n induction matrices -- so it costs ~0.1 ms per evaluation and - builds in tens of milliseconds, while the dense algebra runs in BLAS. - - The 2-D polars are ``ca.interpolant`` tables. Queries are clamped to - the table range first, so the lookup reproduces ``np.interp`` (constant - beyond both ends) and the lift slope is the same central difference of - the same piecewise-linear table that :meth:`_lift_slope_from_ctx` - evaluates. When every panel shares one alpha grid (the normal outcome - of batch polar generation) the lookup is a single 2-D interpolant over - (alpha, panel index) evaluated exactly at the integer panel nodes; - otherwise one 1-D interpolant per panel. ``polar_interpolation = - "bspline"`` swaps the piecewise-linear table for a cubic spline - through the same nodes (degree 1 along the panel-index axis, so panel - rows never mix). - - Outputs, in order: ``h``, ``dh/dv_rel`` (n x 3, from the corner- - averaged polar), ``mu``, ``alpha``, ``umag``, ``gate``, then for - callers that build FORCES from the same lookups: ``h_avg`` (h with the - corner-averaged Cl -- differentiate this for a kink-tolerant - Jacobian), ``cl``, ``cd``, ``cm`` at ``alpha`` and their corner- - averaged counterparts ``cl_avg``, ``cd_avg``, ``cm_avg``. With - "bspline" polars the averaged outputs equal the exact ones. + def _casadi_polar_lookups(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): + """Clamped CasADi lookups ``alpha (n) -> cl, cd, cm (n)`` of the + per-panel polars, shared by :meth:`_build_casadi_newton_function` and + :meth:`_build_casadi_panel_force_function` (see the former for the + interpolation). Returns ``(cl_of, cd_of, cm_of, smooth)``, ``smooth`` + True for "bspline" polars. """ import casadi as ca @@ -1089,16 +1046,6 @@ def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, f"(panel {i} violates this)." ) - v_rel = ca.MX.sym("v_rel", n, 3) - x_airf = ca.MX.sym("x_airf", n, 3) - y_airf = ca.MX.sym("y_airf", n, 3) - z_airf = ca.MX.sym("z_airf", n, 3) - chord = ca.MX.sym("chord", n) - width = ca.MX.sym("width", n) - stall_angles = ca.MX.sym("stall_angles", n) - stall_angles_neg = ca.MX.sym("stall_angles_neg", n) - l_gamma = ca.MX.sym("l_gamma", n) # (L gamma), a parameter here - shared_grid = all( table.shape == alpha_tables[0].shape and np.array_equal(table, alpha_tables[0]) @@ -1154,6 +1101,75 @@ def lookup(alpha_vec): cl_of = make_lookup("cl", cl_tables) cd_of = make_lookup("cd", cd_tables) cm_of = make_lookup("cm", cm_tables) + return cl_of, cd_of, cm_of, smooth + + def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): + """Build the CasADi function of the PER-PANEL section physics that + :meth:`gamma_loop_casadi_newton` assembles its residual and exact + Jacobian from. + + The circulation residual whose root the Newton loop finds is + + R(gamma) = (I - diag(mu(alpha)) L) gamma - G_raw(gamma), + + exactly the fixed point the ``base`` and ``anderson`` loops iterate + towards (``gamma = (I - diag(mu) L)^-1 G_raw(gamma)``, see + :meth:`_regularize_gamma_target`), including the Li/Gaunaa viscosity + ``mu_i = max(0, -k S Cl'_i / dz_i^2)`` gated on any panel being past + its stall onset (``mu = 0`` without artificial viscosity). Everything + that is not the dense induction ``v_ind = AIC gamma`` is local to a + panel: ``G_raw_i``, ``mu_i`` and ``alpha_i`` depend on ``gamma`` only + through the panel's own relative velocity ``v_rel_i``. So the symbolic + function takes ``v_rel`` (and ``L gamma``) as INPUTS and returns + + h_i = G_raw_i + mu_i (L gamma)_i, c_i = dh_i / dv_rel_i (3 values) + + plus ``mu``, ``alpha`` and ``|v_rel x z|``; the Newton loop then forms + ``R = gamma - h`` and, by the chain rule through ``v_rel = va + AIC + gamma`` (with ``dh_i/d(L gamma)_i = mu_i``), + + J = I - diag(mu) L - sum_k diag(c_k) AIC_k + + in NumPy. This keeps the CasADi graph O(n) -- independent of the + n x n induction matrices -- so it costs ~0.1 ms per evaluation and + builds in tens of milliseconds, while the dense algebra runs in BLAS. + + The 2-D polars are ``ca.interpolant`` tables. Queries are clamped to + the table range first, so the lookup reproduces ``np.interp`` (constant + beyond both ends) and the lift slope is the same central difference of + the same piecewise-linear table that :meth:`_lift_slope_from_ctx` + evaluates. When every panel shares one alpha grid (the normal outcome + of batch polar generation) the lookup is a single 2-D interpolant over + (alpha, panel index) evaluated exactly at the integer panel nodes; + otherwise one 1-D interpolant per panel. ``polar_interpolation = + "bspline"`` swaps the piecewise-linear table for a cubic spline + through the same nodes (degree 1 along the panel-index axis, so panel + rows never mix). + + Outputs, in order: ``h``, ``dh/dv_rel`` (n x 3, from the corner- + averaged polar), ``mu``, ``alpha``, ``umag``, ``gate``, then for + callers that build FORCES from the same lookups: ``h_avg`` (h with the + corner-averaged Cl -- differentiate this for a kink-tolerant + Jacobian), ``cl``, ``cd``, ``cm`` at ``alpha`` and their corner- + averaged counterparts ``cl_avg``, ``cd_avg``, ``cm_avg``. With + "bspline" polars the averaged outputs equal the exact ones. + """ + import casadi as ca + + n = self.n_panels + cl_of, cd_of, cm_of, smooth = self._casadi_polar_lookups( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + + v_rel = ca.MX.sym("v_rel", n, 3) + x_airf = ca.MX.sym("x_airf", n, 3) + y_airf = ca.MX.sym("y_airf", n, 3) + z_airf = ca.MX.sym("z_airf", n, 3) + chord = ca.MX.sym("chord", n) + width = ca.MX.sym("width", n) + stall_angles = ca.MX.sym("stall_angles", n) + stall_angles_neg = ca.MX.sym("stall_angles_neg", n) + l_gamma = ca.MX.sym("l_gamma", n) # (L gamma), a parameter here def cross_rows(a, b): return ca.horzcat( @@ -1264,6 +1280,85 @@ def _casadi_newton_function(self): self._casadi_newton_cache[key] = fn return fn + def _build_casadi_panel_force_function( + self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None + ): + """CasADi function of the per-panel section forces, the force law of + ``compute_results`` (``VSM.core.panel_kernels.panel_forces``) on the + polar lookups of the circulation solve. See + :meth:`casadi_panel_force_function`.""" + import casadi as ca + + n = self.n_panels + cl_of, cd_of, cm_of, _ = self._casadi_polar_lookups( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + v_rel = ca.MX.sym("v_rel", n, 3) + v_qc = ca.MX.sym("v_qc", n, 3) + x_airf = ca.MX.sym("x_airf", n, 3) + y_airf = ca.MX.sym("y_airf", n, 3) + z_airf = ca.MX.sym("z_airf", n, 3) + chord = ca.MX.sym("chord", n) + width = ca.MX.sym("width", n) + rho = ca.MX.sym("rho") + + alpha, umag = panel_kernels.relative_flow(v_rel, x_airf, y_airf, z_airf, xp=ca) + alpha_dir = panel_kernels.direction_angle(v_qc, x_airf, y_airf, xp=ca) + cl, cd, cm = cl_of(alpha), cd_of(alpha), cm_of(alpha) + force, moment_local = panel_kernels.panel_forces( + alpha, umag, cl, cd, cm, chord, width, x_airf, y_airf, z_airf, rho, + alpha_dir=alpha_dir, xp=ca, + ) + return ca.Function( + "vsm_panel_forces", + [v_rel, v_qc, x_airf, y_airf, z_airf, chord, width, rho], + [force, moment_local, alpha, alpha_dir, umag, cl, cd, cm], + ["v_rel", "v_qc", "x_airf", "y_airf", "z_airf", "chord", "width", "rho"], + ["force", "moment_local", "alpha", "alpha_dir", "umag", "cl", "cd", "cm"], + ) + + def casadi_panel_force_function(self): + """CasADi function of the per-panel section forces of the current + panels (the ones of the last :meth:`solve`), cached per polar set like + the ``casadi_newton`` residual function. + + Inputs (numeric or symbolic), ``n = n_panels``: + + * ``v_rel`` (n, 3): relative velocity at the 3/4-chord control points, + ``va + AIC gamma`` -- sets ``alpha``, ``umag`` and so the polar + magnitudes ``cl, cd, cm``; + * ``v_qc`` (n, 3): the force-direction flow. For the quarter-chord law + (``is_aoa_corrected=True``) ``va + AIC_qc gamma`` with + ``BodyAerodynamics.compute_quarter_chord_AIC_matrices`` (or + ``panel_kernels.quarter_chord_velocity``); for ``False`` pass + ``v_rel`` again; + * ``x_airf, y_airf, z_airf`` (n, 3), ``chord, width`` (n), ``rho``. + + Outputs: ``force`` (n, 3) [N], ``moment_local`` (n, 3) [Nm] (the + quarter-chord pitching moment along ``z_airf``), ``alpha``, + ``alpha_dir`` [rad], ``umag`` [m/s], ``cl, cd, cm`` (n). + + Same as ``compute_results`` per panel WITHOUT the attached-trailed + vortex force (``panel_kernels.attached_trailed_forces``) and the + viscous drag correction. Taking ``v_qc`` rather than the matrices keeps + the graph O(n); the caller forms the dense product. + """ + import casadi # noqa: F401 (clear error when the optional dep is missing) + + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + key = ("panel_forces",) + self._casadi_newton_key( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + fn = self._casadi_newton_cache.get(key) + if fn is None: + fn = self._build_casadi_panel_force_function( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + if len(self._casadi_newton_cache) >= 32: + self._casadi_newton_cache.pop(next(iter(self._casadi_newton_cache))) + self._casadi_newton_cache[key] = fn + return fn + def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: """Newton solve of the circulation residual with an exact CasADi Jacobian, globalised by pseudo-transient continuation diff --git a/src/VSM/core/panel_kernels.py b/src/VSM/core/panel_kernels.py new file mode 100644 index 0000000..f074f61 --- /dev/null +++ b/src/VSM/core/panel_kernels.py @@ -0,0 +1,236 @@ +"""Section force law of the VSM as ``xp`` kernels. + +Pure functions over a math namespace ``xp``, either ``numpy`` (the default) +or ``casadi``, so the same formula serves ``BodyAerodynamics.compute_results`` +and a CasADi graph (``Solver.casadi_panel_force_function``, or a caller's own +graph such as AWETrim's CasADi trim). The formula lives here only; callers do +not restate it. + +Conventions: + +* Per-panel vectors are rows of an ``(n, 3)`` array; per-panel scalars are + ``(n,)`` arrays in NumPy and ``(n, 1)`` columns in CasADi. +* ``x_airf`` is the chord normal, ``y_airf`` the chord (leading to trailing + edge), ``z_airf`` the span direction of each panel (``Panel`` axes). +* Angles are in radians. ``alpha = atan2(v . x_airf, v . y_airf)``. +* ``umag = |v_rel x z_airf|`` is the speed the 2-D section sees; the dynamic + pressure of the section force is ``0.5 rho umag^2``, not the freestream one. + +The force law (Gaunaa, Li & Pirrung, TORQUE 2026): the MAGNITUDES of lift, +drag and moment come from the polar at the 3/4-chord control-point angle +``alpha``; the DIRECTIONS come from the direction angle ``alpha_dir``. With +``is_aoa_corrected=True`` (TAT3, VSM model) ``alpha_dir`` is the angle of the +flow at the quarter chord, ``va + AIC_qc gamma`` with the matrices of +``BodyAerodynamics.compute_quarter_chord_AIC_matrices``; with ``False`` it is +``alpha`` itself. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + + +# --------------------------------------------------------------------- xp rows +def _is_casadi(xp: Any) -> bool: + return hasattr(xp, "horzcat") + + +def hstack_columns(xp: Any, *cols): + """``(n, k)`` from ``k`` per-row columns.""" + if _is_casadi(xp): + return xp.horzcat(*cols) + return np.column_stack(cols) + + +def cross_rows(a, b, xp=np): + """Row-wise cross product of two ``(n, 3)`` arrays.""" + return hstack_columns( + xp, + a[:, 1] * b[:, 2] - a[:, 2] * b[:, 1], + a[:, 2] * b[:, 0] - a[:, 0] * b[:, 2], + a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0], + ) + + +def dot_rows(a, b, xp=np): + """Row-wise dot product of two ``(n, 3)`` arrays -> per-row scalar.""" + if _is_casadi(xp): + return xp.sum2(a * b) + return np.sum(a * b, axis=1) + + +def norm_rows(a, xp=np): + return xp.sqrt(dot_rows(a, a, xp)) + + +def scale_rows(s, a, xp=np): + """Per-row scalar ``s`` times the ``(n, 3)`` rows of ``a``.""" + if _is_casadi(xp): + return xp.repmat(s, 1, a.shape[1]) * a + return np.asarray(s)[:, None] * a + + +def unit_rows(a, xp=np): + """Rows of ``a`` divided by their norms.""" + norms = norm_rows(a, xp) + if _is_casadi(xp): + return a / xp.repmat(norms, 1, a.shape[1]) + return a / norms[:, None] + + +def _atan2(y, x, xp): + return xp.atan2(y, x) if _is_casadi(xp) else np.arctan2(y, x) + + +def _matvec(matrix, vector, xp): + if _is_casadi(xp): + return xp.mtimes(matrix, vector) + return matrix @ vector + + +# ------------------------------------------------------------------ section flow +def relative_flow(v_rel, x_airf, y_airf, z_airf, xp=np): + """Angle of attack and section speed of the relative velocity ``v_rel`` + (``Solver.compute_aerodynamic_quantities``): + ``alpha = atan2(v_rel . x_airf, v_rel . y_airf)``, ``umag = |v_rel x z_airf|``.""" + alpha = direction_angle(v_rel, x_airf, y_airf, xp) + umag = norm_rows(cross_rows(v_rel, z_airf, xp), xp) + return alpha, umag + + +def quarter_chord_velocity(va, gamma, AIC_x, AIC_y, AIC_z, xp=np): + """Relative flow at the panel quarter chords, ``va + AIC_qc gamma``. + + ``AIC_*`` are the ``(n, n)`` quarter-chord induction matrices of + ``BodyAerodynamics.compute_quarter_chord_AIC_matrices`` (trailed system on + the quarter chord plus the curved bound-vortex correction), ``va`` the + ``(n, 3)`` panel inflow and ``gamma`` the ``(n,)`` circulation. + """ + induced = hstack_columns( + xp, + _matvec(AIC_x, gamma, xp), + _matvec(AIC_y, gamma, xp), + _matvec(AIC_z, gamma, xp), + ) + return va + induced + + +def direction_angle(v, x_airf, y_airf, xp=np): + """Angle of the flow ``v`` (``(n, 3)``) in each panel's section plane, + ``atan2(v . x_airf, v . y_airf)``. Applied to the quarter-chord flow it is + the TAT3 direction angle of the force law.""" + return _atan2(dot_rows(x_airf, v, xp), dot_rows(y_airf, v, xp), xp) + + +# ------------------------------------------------------------------- force law +def section_force_directions(alpha_dir, x_airf, y_airf, z_airf, xp=np): + """Unit flow and lift directions of each section at direction angle + ``alpha_dir``: ``dir_flow = unit(cos(a) y_airf + sin(a) x_airf)`` (in the + plane normal to ``z_airf``; drag acts along it, TAT3's local inner flow, + not the global span axis) and ``dir_lift = unit(dir_flow x z_airf)``.""" + dir_flow = scale_rows(xp.cos(alpha_dir), y_airf, xp) + scale_rows( + xp.sin(alpha_dir), x_airf, xp + ) + dir_flow = unit_rows(dir_flow, xp) + dir_lift = unit_rows(cross_rows(dir_flow, z_airf, xp), xp) + return dir_flow, dir_lift + + +def section_force_magnitudes(umag, cl, cd, cm, chord, rho): + """2-D lift, drag [N/m] and pitching moment [Nm/m] about the quarter + chord, ``c_x 0.5 rho umag^2 chord`` (``chord^2`` for the moment). The + coefficients are the polar values at the 3/4-chord angle.""" + lift = cl * 0.5 * rho * umag**2 * chord + drag = cd * 0.5 * rho * umag**2 * chord + moment = cm * 0.5 * rho * umag**2 * chord**2 + return lift, drag, moment + + +def section_force_vectors( + alpha, umag, cl, cd, cm, chord, x_airf, y_airf, z_airf, rho, alpha_dir=None, xp=np +): + """2-D (per unit span) lift and drag vectors and the flow direction. + + ``alpha_dir`` None uses ``alpha`` (the ``is_aoa_corrected=False`` law). + Returns ``(lift_2d (n, 3), drag_2d (n, 3), moment (n,), dir_flow (n, 3))``. + """ + if alpha_dir is None: + alpha_dir = alpha + lift, drag, moment = section_force_magnitudes(umag, cl, cd, cm, chord, rho) + dir_flow, dir_lift = section_force_directions(alpha_dir, x_airf, y_airf, z_airf, xp) + return scale_rows(lift, dir_lift, xp), scale_rows(drag, dir_flow, xp), moment, dir_flow + + +def panel_forces( + alpha, + umag, + cl, + cd, + cm, + chord, + width, + x_airf, + y_airf, + z_airf, + rho, + alpha_dir=None, + xp=np, +): + """Per-panel 3-D force [N] and local pitching moment [Nm]. + + ``force = width (lift_2d + drag_2d)`` (see :func:`section_force_vectors`); + ``moment_local = width (cm q chord^2) z_airf``, the quarter-chord moment as + a vector along the span axis. ``alpha_dir`` None gives the 3/4-chord + directions (``is_aoa_corrected=False``); pass :func:`direction_angle` of + :func:`quarter_chord_velocity` for the quarter-chord law (``True``). + + Not included: the attached-trailed vortex force + (:func:`attached_trailed_forces`) and the viscous drag correction + (``BodyAerodynamics.viscous_drag_correction``, NumPy only). + + Returns ``(force (n, 3), moment_local (n, 3))``. + """ + lift_2d, drag_2d, moment, _ = section_force_vectors( + alpha, umag, cl, cd, cm, chord, x_airf, y_airf, z_airf, rho, alpha_dir, xp + ) + force = scale_rows(width, lift_2d + drag_2d, xp) + return force, panel_moment_local(moment, width, z_airf, xp) + + +def panel_moment_local(moment, width, z_airf, xp=np): + """3-D local pitching moment [Nm] of each panel: the 2-D quarter-chord + moment [Nm/m] (:func:`section_force_magnitudes`) as a vector along the + span axis ``z_airf``, times the panel width. ``(n, 3)``.""" + return scale_rows(width, scale_rows(moment, z_airf, xp), xp) + + +# ---------------------------------------------------- attached trailed vortices +def attached_trailed_forces(gamma, va_panels, v_induced, legs, rho, xp=np): + """Kutta-Joukowski force on the attached trailed (AT) legs (Gaunaa, Li & + Pirrung 2026, Sec. 3), ``rho (gamma_j - gamma_{j-1}) (v_rel_j x l_j)`` on + boundary ``j = 0..n`` with ``gamma_{-1} = gamma_n = 0``. + + ``v_rel_j`` is the mean of the neighbouring panels' inflow ``va_panels`` + (the tip boundaries take their one panel's) plus ``v_induced``, the full + vortex system's induction at the AT force points, ``(n+1, 3)``. ``legs`` + are the ``(n+1, 3)`` leg vectors, bound vortex to trailing edge. Returns + the ``(n+1, 3)`` boundary forces. + """ + n = va_panels.shape[0] + to_boundary = np.zeros((n + 1, n)) # gamma_net = to_boundary @ gamma + to_boundary[np.arange(n), np.arange(n)] += 1.0 + to_boundary[np.arange(1, n + 1), np.arange(n)] -= 1.0 + # boundary inflow = mean of the neighbouring panels (halving is exact) + mean = np.abs(to_boundary) + mean[1:n] *= 0.5 + if _is_casadi(xp): + gamma_net = xp.mtimes(xp.DM(to_boundary), gamma) + va_boundary = xp.mtimes(xp.DM(mean), va_panels) + if isinstance(legs, np.ndarray): + legs = xp.DM(legs) + else: + gamma_net = to_boundary @ np.asarray(gamma, dtype=float).ravel() + va_boundary = mean @ va_panels + return scale_rows(rho * gamma_net, cross_rows(va_boundary + v_induced, legs, xp), xp) diff --git a/tests/BodyAerodynamics/data/panel_kernel_reference.npz b/tests/BodyAerodynamics/data/panel_kernel_reference.npz new file mode 100644 index 0000000000000000000000000000000000000000..897a218e9f7388822722217b06d42415f251e143 GIT binary patch literal 215246 zcmeF42|QI%*YM31Qc;EmA<;n5K!|-pQihDBBAF>uNRy&ah)fX`A~edB3}GKCG-xuU zB2*eBnHq@i=A3hL4|?D4?ft&z-hTCae$RfaXP>psT6^u`?6d#(=xMRCNi$*pM495% zF4hd-|0Xj{W^(s%vv9Dp_wX`u+hOCSv1@^|>uzSIiA(|b6+{`SYpkQUftksZX|JTU zjk}ebq@sxAa@%E+QX-PJE^ZVp&gL#|);2%o8WxW3HWazLy@jg{#s<@f&aGHq&Y!ew!Y^v76LyVs`lU{l*)*%`~jvj3u>= zufXgu!W;f!W16Fbtp~mKoUsVIcrPILoSlnYxJUQL+Reg_Uj0$;s=IxED?%CxzMqRV zB>a&W);>S|F?Vpb-eF~9ZSLV_V{7X{xA`aAXz*|#WiH9Tp|jBN{Bw7(wxQSd8hIX4 zKO_!__8PuMuyUkVZzOisnk=?_3TQ{}?9j8ze&dZby?V2K+9I2SlHnw8ymp&u|Hd14 zdiAC}Zg%VIToKZW$C=YBf$(O^u#v&aiEbAOZ#LRopHvl5NSkWw6~4^)ojY{9N>Ccd z&MYNy8`8StFlp8O-?>A#!}LYRWtLC*Hjr9J*e-V>;g0C=WOFpPcAz{BIau!SaBy)Z z`VkW|HbxTVA3X=j%1QmP{YzZ>J8wUbV$sElm^h_;-A#;umOC;OF~Z6L-Ffua6N=x!VY1VHbcC? zYWSOwS~EmQ#}}2#xI796t_2Tg3)n<}y4kXwkKTj>#joxH+w{YMcrVX`$JJrLTk&D_ zyggw+-&wh<1%!bQ@1ngA&k6%?@lp$y`{DIw;nyL!9zw4RZilc-8`8OJgNcX^jF6@7 z3=vYi&-}*<+K}$qEzd{9RAgAk*BCjXEcngi<)8L@k;Q-1^En0YHvM?k`Ua-dBUiLl zCLnk7BI_bU5hO?JgvMpL@6gMzRZ>>D5E!33SL5Bk5|RDHzh~(lW^|sEolesBaM@($tijSE&LfdkL?g>hcnMU^)%`Fgc7(+{*Q4rT7T zU_q6au3szn%-Ip>1}N3S1TCt(^U2IHxi!YXEKc;>iK|Pg@`dVt0xe-Hz&fsFj%yW! zsdC<~pQ$t#R4!T5c5TUb72-Zu6)yI;m%Xvj&X)48OF@B<&;K%5AI@7$TVGdoP(SP# zN2`BN?cn{uJktHl!|oq$BHe!%<}Vzxgx22yEI!bO6rVF#JYgOwo|;(v(KZWQ8;q3GJ329qLZAM;Z?B=RgaD>i5C(?(v6C_&!?LqHyypF>^z}> zPO37Qb?6yCGDXF;!nA8Wve~m^jMg+o)ZV3FPpcIl(j3CYuD)L%ISR!d@orN@8^#=- zyhnCC(r1%za#M9B;>$eF{AR2onl51WXg&`MqTZjI6YVR72)+^evNc^1t?9gSJ<8!N zSdPij$+-LuCWp8ll+bep(}Uu6QNnH=W*6nf?h`e|@6!gmPc#qn3%!N=C5ZWj;&DKU zIA~yTK=HVtMBKcwxS@EQQBs^qn}Kqgo7{yQ#SgoY^aZz zK19dD%Rd^wWd3wzEW%0KPlkFW`CG4?9IRbk9GpEE#4+f+^kTz(6$cL~34Bd%vPuJ<)=hp?;?jxS>YI+bsKP{c{#Z!3_@d-4?c+Lkf>qe>j` z%`Ikn8XpHFc+Bbz?!d7*k1@KXm4P z%vHHH0B*Pz;3F*x%}=Nf!|sNwu`+u z0=w~&kQ4Pd@aybxJ%nDu4q=xzPc62V{ch5CVvNDm!DGF1cefqr=FHoU#@d2iVWx2j z3miZePk61$3p-%k&nIef$o|=$koVQ~ z)%F>e!9?XBjyl&f!0gXS-#rU50sqd2fz0MvAcu3@Jpa2{U|ghh!8^Wez~#NNVFJho zZYw-nw;BE_d2sm{ydH7gWL%E~u6Gh{hptS@ID5n(4R{%R@J%|H0qSq({-`;Y31*+XSSX>K z1vrk*IpL9#1=600UJZJY1=3yS?s(px1uozvAt&n1!mn$_^$>aqJA_@@p#FY07|T&R z8y6=V%8QHtXLxVBt&?k^)Y)LUPSeIsZmfMFEW21 zSJN?t7nz5bgj@u#CyHP91lL38CG2G5c4?y;Y5d5Ni0be+PZn-AliR!yE)~cD807N;8y!-XVM(h?Dk^<;gH}khbSZU5RYQ;k&r0=Xzp43=dV^q8;EI} zj5~10VZy`V%QAe)>0(RiU^edVJKU(JVhJya@;~h@TUpFn$6bhoy9!MzWYAM<$ZX`#N=@=lz-HL{?``a?SR&`tJ)digljYt z)F3ZL$Lv{x9cX(re+5eG&d^i;Be_Gj^DUOZWcV!Sid4$)T8C~A`;PHq>yaA9(5zga z;PM-H=p8@PCLvDqo5iTl<|8DV4&5f(y1jLRy~Fr3rt+(Ye%Cs5n{0!9$7d3$b>0>! zOeWl+S*v0r8$AZ8Mm?L)TbWG$^NlfBG_zx! zxl3dpf@yJ^?M`>rg8Glq{prn5K$q*>w)w}Nfn$^UPT!5M2eqr$$9z+30CG)hveJ_q zfUL^Lyh9WgZ7yCCa-v>1eqA4~htNydA?(sVk5!;Q*8DD}BSRR(u>+#~qw(o~f$x)? ztFU|3u4u4q%r%|~4pD$3DMC}-dpX6 zJf+#;!TLU9;AnXt24uN5R8Vm==d8m9$HNfBeLScvaOtArXhIsPEMu88kOPLv73UWW zKmtb_|Il9jUAZPAKYr$qRRb(k98Ev1NtAaFDy4n*7t9F2$@Db6Sn zXWEn+ilZ$W`f7nlsee2g6O+$B9u3>Q_}}HxQsy2F`OR>~7(iOY*g%okiRQ@X;Hhq$ z$xG0J({CkWW{pPzlO(mLG+7`Oaa9K2DL~qr0mY5ppC%$_7dx7GjI%)Etn=^BxU~$u z&2h2w{DCP*AA8nDG{X!LNf2MwNdeMCZk_bEvz~%zCgkRDbgxGeracn9L;=#|@|>W0 z2rtqsCn0|JCy>UsOxThFq_MnzVDkxZBMGaoZ;(=3iJWY+=Hxt~h)#|PtWk($M(VYr zB>4TL5bL7>=@%$Knn}s7gDsEVfND&RZo%amm>lAIP?8=Lw~LbOqARicMC0-Ml*R57 zW#<_53l+!xT8#OH;&DJpaX|67p`^H>c$`rp&eyOwlXljJqFb|v4`+V^kw$X_0UM=< zy2W5fT9ahTH)0H+Ir@t~&M*j*mThV`uEkXdu6T9$e8R?+&@;R2M!*_9xS!)lN8%Gh zsABMAhT`mXFkZznV%n+ou%^>``Qg1LFl@&=$*HeQ;C{R$#}e?gkHi9 zVV5>jhf--Yrz`R9{|`(1g}43pu(UM3#!J5&P#J_v+d1vvxhlmBuyDeSeM00OQI98$)c2hjV?ce-q7$YW?hhs8XSc_03xHaFM=J*? zV44ZkWL;4b1ht#5Y%duE(^%rUP9A~0H|O2GLjlu{q*ig*P{6dqcuB~&&By9#;n$hs zdI-IQ9l|bc-MhT9*@8s<^kmJA zwjiQZL)FL93Y=QrZl2_429__Xuv9WL0B^gC_vF7?25PVp;&O=BgTdH!FbC5EalH_? z0~0X2q)Z)(QV9=V05brZM$k+I-9*q%{~3yEXgq0->R}TU7SDfd8u*{HX(Mb7|AW(h zzblOYzqG^YP&jSwXg!~?;H}@$X(K%b|6O?62>TcR__ROg(?%L^{uV%OgqK28B!$5V zwGo=_AENdTb)!s${}u8_J53Iyxuy@FCjSCPZKTQnFN4%ZILO1$lo}1xg*{m?IHo2v zyledpoZ2Ph2a&|?+6e12#sSrKXP%!@me&YQb!UnARAi0iRKkp;hbkx8rOI#6-X~S=OzRg)Jb#knK$Vl?MwOG| zOqJ8l!9z*45!Q7KjjFw}^r4>Cp7e7N+-c)V_n$~R z@YmM18=hv>q(kS4580@{<4$&kZVv!^txjEnmhLCs-PFNEZZ>Y|NQ7Kq#`bc>yR-LI~?^tu!&!J0FosV5S6X$NUKM$QUE zrdB#R#glxIAA7i%$gp9)(Cs=vZrk$K9mn31e9^n9)HJFubh{7;NLq$&Z8afvmmlFz zso$gZ9XaH9*y*RU3W$6tGE|_askl0Nkgqq)z+V0EExJNT`&3 z4vq?Ggf8iA03J5QJm|&-z&qu-AFE|bU66d_{lZ`$KonY zH>*aV#*{fQZAK$#G2@5+V$T8c{rT4qMbAK*S@SOIi#4QpTEmx|L~fAd$ry+k^Tx{u zS@ivY;LOl>+rt7torj6$s#}2|H>}cv%QXnBGw=4Ex#kGSGn5yeR~ihagx~#8j)sEA zyVf5+%M%S=Su)QHKX@FBtKsy2E)Wjl^|rXS@q~iQZD)4wO6k_Mfn1nfAqAew_U%0apn*9^(;k@IV&u#CuXyv_vgz<>lXEZ z6kWQ}0ilZ>ydsN}Nq@xl|oK; zm0Kn9(QwKa7&D85+p_>c-Z?ocIcTz?LetGAj@#Y?TPuRmQ%VT3dG?(ddcMr)2?>eA zXMMXN-#)AJ(OfD>U0g=<+cV!G*LCU19QP-o@!`^q+t6jC1^vZ|52%mC3uu=p8kAp0 zXWY-6F&Ttj)vQ~_jTl)@S2OAUOx;B-@YS0wx~UY1DTpnU)vcz=-+q`92mK;kwXy&jhf-K)%4I-YKNG96%+2!(dVt08@S+oiiB$PD#cZx|MzY{* zxOm6S-Ts8!XSe>=Wp)X`pN`D|Wn_Erd)^G?n+`R6{0Q%S%`hx&k%b$p1~ab`sKYXl#?d_fn)$ zcl*pW)(WV==lW^hpV*METMyk8&zgd`m=ql{iCKt#J74S5dEh-b_2hx?13`Xt`er_(zf`pzMV4+XXdI)Des%NgYvK|| zkHvBO^ro*vCHiG|8!yyDBTi~}1Su>g?d%Oj$_CHZ{F)Vqa#xJ{fwRFwts}d{U?|y` z*^Y`tzaRX<(J!wU<y)!sXS zPc31eX0eoqwk6Ch;Hj!&vVh!A?0!6`+yL8JZa+S!wgINA+rPe@zXg_Xd{I|nwS*R$ z`KmF&mhio?;E%2#OL+7svibR23)r2mbaI^f7MNB3T4Q;!328m}@&1l4qXv{Qh=*pu za{GC_+X32~>-65!2xxnlimP<_tywd!jOVNxX9Tlw2Yuz^49^!aVZFywi0ueYl)WJh1= zzI|*2-zECLUWE30|6)RP~0P055h} zcQ}4H1@>OcR&;n42R@iz+wpXIJTS0dGb`j{9B|iU7A{;D2lO{suVz;{1p?&H&dD~R z{M~9-C$l3SbnEgkMNUfqFXzroYt>HxKlXNbe&yP?I?0RjrLv0-(do8)1yFczTAKNc0+_U4JZpm9aP%BgbXOwj6-kMwh*Ir7q@umDd|622Q zzUc+9;p$zhT1vSPyqd=s<@PmFM!+D1UaR5<&)C0&)k5Dj*-)o zL9ns0X%m~1ZeIbtd$;ZBHp&ApIs?5VH(v$e6;Zh-E3*LqvXA>Gf;8|zQcvpUwq)?> z;>UyHJjX#|PKZ&e`c7clu>L`5)=qc^-NoZldmQ%qJeqV_It2>esIgdYei_EzXnDCL zI2-C-)atL-$%RwTnH;(rk`FER`f9x}$|D6-=dmsIV+q;iSOb95qQWOjPWesnSV!RAkv<*V8lM0ngR# z)oNMrtc{j+?50doK;;z-P7wWlSU?$L0R7B%6!Pgrft$M(=kaqM0^O7oqxfucy@Kufh24B4v}-8ID) zE)bqB;8pDg-^#q@ni9ANdVL4+ERTKQ67~rqIiNB=F|R%wp{~Ds-K=#khun2 zD>hnA$t(c#E!#~K;_`u@)5Xx_=C1pDVt^wLpO<%tYB@@Xx$TI($j2H$J z%e*D%qg0RuLT>NA#>9MZSy*9AuU$UK z7I}N7ATAF)lol;0$+!yI16XX-ZfBArI&W&qXfGfb1QEkNDx`v7d7+H-5B`tgCF|f; ziwZ(6iTZ?`_?sxPPc#40PuL;s(q`zL^y62F3>`eZO3qM5!NVNY+X`}71xXu@o6fwi zBHSZC)-n(@14S%xsfLin9L(ib0b&IrKjj$<3ikKUa8MMG7fn$&A=l3 z2P-LIOtv~4Boc=70*OJsS@vV__Ggm(8p^nkUgR(`4%X>hK>d;}sks;g1v^%>*GuD} zG-9*aMA;o}0^Q*O;fnmzDWeWc9YEXSPiWisuP7hR)YtGn|mGsg+E*4}}M$f7bw&SRBb0n^J`?>)jvEY^!&vtgIe9Gxh zX^9uV1Fw>I6-ge2JTHfjfH(RLkSt?^5CdM^^5_W6w{YDv8eskQsof@T0APEr)-{@KfiTjrg zT545Z_YMe~NYu+^wNT~SPYxJ#9QX|SyWOpqJ@25(Ju9+>HwpIuuRS?YTN=MpjD&Z+J!4(wFrH0~Nq;(oYqh0gV$57Aj{hEF-TvuZky?#kmmC|_7rM={m_@y(%c8%Lc9NpI zGN>EXzm}fZI7ZKems0ba(?%yrM&Ey*6F67CR1MN{3qE) z^@VORFSS?|pVo(Ar14F!We(G*zR)e_^~3tiw=I6Mq#Rq56;Ay$(Jv|z8xQCh^h$BH z=h$IS2EByJR&5~~X1FtSyCUehq3&7{LnghtrY*#F-zgpqcZP141h?D5o*&r0iPU9+ zWV#*^E%WPO1)gS%!Fp}ja#2usObxJlt_NS&J_VEWZM1e$STECKkzErgtXCWocdnem zdR@E!^otCI_3~F-BTz+Qy-c>d7QRxd2Yq}Gd!nrxz=9hUPlb;@2cArx&3h@VS80n= z`UVQ?)&JN&_`~!@;Cg4f#3TyqrQhyd?pyc_@Lz9W9g|9aCh@AP8{w!Q9#F=By%z4N zITBuS48*sqb#F-w2jVK;+tw@z16}oI(0tfPup$(2@m*vI ze;)>r%U|X#=7|L7k^Sq|Eh8WG`#yF~?x-2c7_?U+Z{nUqTfBkc+qj5%oS|Sr``iT= ze?)-2wkF2|-H(A*L_m4t z2bX+S9)sH>!_M_Di-NtmuZ}J7je=bF9y=}NkAlHZzfMexi-a~zVf)U@g~2gf1Y=@@ z4v=OR_9k8>BA~x^X|N1s4BSiR#&`Zl?Vq4`Z=y(d_C$nobn>a3#ZqM3oM7P2+Yi6- z!Hq6wJHS1evs1U6S&oc5+!kTBgayr%HZ2d9=>`S|%MBYP0a6uFQPv*7issZ*y-SGd z0S`>-Bgaifkj(FC0x~aH(TR<13*U720Qt+s>r&+qq}BJqeVuStH2Z}^{sl_?CogRr zaw)jij2q?(EXgeB4b6(|e(!E@=S20Vh$w)BCZ7@wVrNG8`ETakd8QLoArHScnX4eq zd$^aF6*8ee4vGvE>2|`PT^kc}DE(k+&(oRx8(Gi}X7DIvem8U}3vghj;9gPrNhKT9 zSkdUUj~mV_^}sV7JfdqTxL4(j?z&k)tSEcbS(lw*J@D>z*Ch4^1YxiFrqwUXih5+# zymX)012@GVnOS-bK@x3rL1Y0l+Hyq4{cK7XyuGWEw^tcK7P+hvUiiHaIz~KSexAJp zE~?B_Q<QAGz--q&N2N z8q{>HV3=@$FtS+FKgUEv3u!QO`x<4eflfN}bt7+*Ao3{s@Q+&*;42ZlsNYEez64F0 zJ;&uvLfTWOX9wm~kRLwPfr zCq`0$FS)1*yh4Vms78=ii-+qZ)VR>hV5Xrqx^&OjBZ7LX(8jWL;!~FjqMyZ5p7<#l zpg$s%^V|EhP&I_V<&2;(y1%vl=8qy{R8t7t7ie3H=C5~pT>DuBHM8F1aOLTGG}>`d z;BIbxw6wbI#^UoL=mL@P z2IA}2!GPgDF&Of-{(&Yrw==OeD(!${a=5M_%zsXP^M8wJzsa!S$8nfD_zWQbalqq@M|T}( zw^9C{F1S^1>{mmu`OX1(S7ifG8qRkWw&}xrw+`2zIBy8+Qp;zqOErS$??_w_>NbK# zJk}m!!A7vEB_l>zn(}eCq)79rMS9RJLQm;Mn-eLbulG!$QZ>VKlrf+M>lGdG+n+^( zvkTODR>htGLZazmO`lEyt|;kKH4SmVLE!B)QQ>${5$s*nMWMcStyWXGVH^kU9h$2u zB6kXeGQU0!^y5JGLj!JQk9csy*&(unLVfiZ^X->5OaN;8+x4xcB!DRK`K#dkcrYVZ z*}du0NzmK>`jIH-anfeRQ0nW~@kYGY|HHoiZ?UiK4F-avc6kQjzBu>GzD#&>1*q1= z>?tzM12@j{FUYXC3MSqPy|CbZ7I@Qo>F`D*4Ky@=W}4`j43??uRqKpD4lazZPYF=k z3HTN!&Yp8&C;V8$-g>g;IOGV8%^|gKvgpL+ymcQ#a}7LJ`Bx zIV+;_;U@3Lv)@{hxvy2}2PTXjQN|cVWv!m4g=wZi?s$t@p|zRtp53^yvo~eH3eEg> zsiI5pW$=@nPZaRWprj)ip@3htMhXhCiYMWNx6CWgeGG$ImCOz9`Qgwz{_Omckj>X!{n4BL}s@fsAsFRrsNUpz(U)cK&k*fyJtYlKgvoKv&&VXvMwGq(9pY-ZF{qUfi=A z?mVm3x__1r99SKux>)KU3z^aUyn4d$;y zbl8makNx^%_KW;R;vf7)nO&5T|2h06F=B`Rd+4uWkGOvq{zcxE_{V?!IsY|mjA4v# zg8kP382R1K-$H7cgMM=J~%20wbThVQ32M z7dM4HI50Q{Hu#OUA>HY3_!lnWFRED9cFoo4ro1oL@lcW*X@1QRaTwh$jef}M-`a9|IG1WWKj1uhRE!Im}% z^#@W&FvUAIM;{I$!Gem!yI3eBnDBa=WsyTjunNv@DW1G{AZhn_w`oI2ure-oNi@C* zIBm-8wJ!Zdf-Pi~t~0#e1o|I&dbT?KBEdFfa~Z#(kYFk9m$r%zA;B&xcx$&(NU(XS z*Pbg4A;BbPJ`r9;A;FNH=R_TckYIW}#wDy266`Z~z*(swB-lswk`Q(Z2_|**|AhUB+*$WP#&f& z+^oEe+;-S_(Ju!kd{+9jzK0y?Ww%mKcJs)N(Cu;D!80W2+#=eQoO8>IZzDTGw?}ae z0|V8awmOjr5SgR~>Kj;^)S_SfYp_Lp#7Hd?Ya&8@uSDwz-KN%se954TnYN@p;+^6; zdBoJxq1)6NY!N1DQj1hA&u}Ch`T6}D>ZnD})K_U=sh_w%xtB1Tq)>nBiN+OrM!wD> z7R5Km%qDes!25>k%FmYv)K+06-s{_97bETW{3O@hbhgWkm~iwAdtvj72&wNpL?X5k zuF#BT^hJN-OP6uRQXdINCJr3{=os))&hYe|AA6kS$(v0RZja&#-Ol-M*!M0h=%^>X9_L89 znlp+gbUW#%WbZoub$=46P2M%X>LWb)`5W2P>4ibQ;T6>9`LIu+98{*iGI*t41(N;E z_}SZ@fUCBB>Kcz4!Nu>LCU)CifHEzI`|ftHfPbO?+y~_?U`GS9y|&kDYUzW^ujBPn z@#{)(J@L5So4B25+%9cj-OROx`YVjI(PRv$i?!R@W$&>t(6(iUzs&4Nus}huRJr*W zC_KR%o^22fR$W(fi4Hjqifgu(z7324r9ReWR~2JHY~SN5-so76h?l2vxfWipvlP4T z9j*u9dhg{%&hcqHj_#)WxPD(;j9Arr+u!eOsWPS-Q&5x zFjWE7{n00qEFORfORI`bi{1i`7w>r(9fd%AYv+%u11UratuQ$h!s^k@f(*sdV zFU0M@49qTR%tGxiP*0lvJkmq$>5PGM9WmeV2CRDs*WOzqdnH~J=_$MLR(u6NA|8`@ zl<7wcoUT|h#}bGkvc2Cwtk)GoOvdd!n8x}U3dfaS-%(BZoUFGmx_ZeR~}1HmFI_CGldd<#;KJub?o z%8k&L{R&BkfWRs<*$XRv-A_TJ>V~jMmLTn$8Lv#5%&2W{gRf2cD6tN&B!P2 zGxFGAefWj8zK;~9A2w^y>c{T~9wXh)cr zb9YRWKvzj1pIic`A{`&&o&>ZV(mPM0JtJ^LuK?a zX>`f7)TYULyhyds&g>IjdtvI~0CTPWg{ZA{=X`dZF-T>p#F|~57BFY;80Ucm85CKv zer)ZxPr&QOpd4Jmk#r5+Q-G~k`7{c|a()ek)n?&{DG*8nFIwh3HCdQ#e;X+j`s; z8`u_LQ@bbF5#Ek7KIAB~9bQ>DcUIv(cbM^L#`(4>J7An{z2~jhJK#yY%);e$c)fc3 zy4ScKQC#m^+)g=emv#W{n08uYlyEWzcd+1{pvbM=(O~iO$3KcSj)UHNjc?zu#ell% z%_|=TM1wTZkIStLBEi;!AtK&KLcsI|yXHvx9t6tye(yS#IuIok#N^Nis|U|w*TLtQ z9>{ZlP%o6h?7)SXUDA^J`?rL9qlA+&po0awSzpLl1%slC<yt)tDPZdR9^p;(XTc8U&w2N^o(Iorr4-8KE>KG+TrP~)JAq%fku00J46b z8GmX@0o3=X*54vh0HxRLy|jEz0aSWFBh_BI0P>&dO8dU>S4qf;dc<{XxSlh(-p#lj z!Y*wrcfGPpBVtKAa54z9<=T2xaH+6x4C)16IS1_k3Cu2OH(9>ERBV)RGRDBEcbTdjVDENhqIN7o12Z5hYj^@7!$_XFcIY+{f@h=oRo;{ zAN`m+=fi%`Mod`?5w;ME#;mA!g!TmL@%R!#VmJ!?|p*;P~bH=3Cc*=(Ta#Kl8ZRSA$;S^0%F$j`%*R*)9)V)~< z(X$iSWeX3h_ke%$|37h{?fCm%1-~M)=_f-%U7*oDkF{#V-XZkfk@{Aac1mX`7)y0T z=Mr8CN)U}jb99N0Nv!26Y;$*>6(@a$Abs_un$K0uDN7oqE;G@ZDUDANusViDZ zw8~FJj_M12GgzaIx3uHK=T>-uaE@lBH9CM5gOBV)`A1K1i>3bP@l8ZTgi~bQBHQ~N zbzr;Cqz%<^4+=4NCXiZjPV%j`5glb58mmn&OF8SrOK|hS)TUYNVcwy)ULHrtVnTo=$)m@$04`RxgRxA zc?Nt#7vAdJ#)McV?0xy%h6}lRYpQwj>Zf4Q$73HfxW5A_o-5Zc%o~sNbMCxvmG%VW z%~V@g($oPwXI$9L@5h7iD@M02So;LH9YpPSHNFPSfgY}G%T=G~Cd3qk0krhXVlC3F ze8Kv92};{OLr=?-JfYjuvx6fmNQ~ZevrC6=vuj%NM^R25GLQ9}Hlf?>8tf-Nw{X_Y})3A=EF`(E38pNET~Xi;wrLIi#i6P$Y}ksl?ia@tLelVN>dvYua|A ze<*7fX-unjy-bCEX~#7BrLq=H3jDV9A(_e=XB*>9j07Z4s4t?34+#_HA3etW=|c?V zF~(?ZiNo&O@uU~&AJ1=-9mNy6<+7wE$T@A2jU|yY)`xCU*?596#YO=-2D1cYcUHt3 zM3Y>RPJ81;Gz!@jx}E*s?e5NC5+jpYziSn`o&HB|6&rIQccN8@BRqJ5F~&EyF2pQa zoIJZ6SRSiYpQTa-#7tWGwBJ1e9jDbll-z3sYxqKOF@= zNp5L8zxX&X)aoqYTN4AGSS-_dY8wkYCg&Bt+j9a;Q42jkV0RMu%*owmH2D)5J3rNZ3wIT` zzlOi<0ADp=c0nWiSgOF$)g1GTtsYRm2hgSL(S1;m7V67M;Qk2H> zMQ#qVGjFH{2eYYiTUCs-U%2|Nzrc-)?1aYN7IakIwb9Dv0c6~yCA+c~EW!Ly|C zEq)@K{_$CqZ%qB;v#>3P|6M+dFWKnf$b%ST0a~kDoyAmcbE7;pCThNc0*Lie4A3&G z)4w0&BZ#s(XhrL}2_hAfFhFbC+O?HIUq#Sw4XHA%_EV828w}7|8aLm0?Sa{-?|YLm zXCDY59iK5kOFUZ3B=h)O)c^Xf=!_y491w-mbI%>EnUA15Kzq8Okx^W8w1 z&!&aw#XZkA2i9^>0WBX*<%>sjWYEI+wHr2U`UPka@*+%51+;i^JxQ1zDxk%>53?(e z*+qkzF+giNexE+reNq9fmAGG#m|s*tYbG9tvCjwNKn1jP@VNEa55^6hivwEg@HkJz z;!N7`9tvpv{Hnvx_mtQm_BRk(#LNES9x)iAHAAdfmFV=uky*?sMq#uvJVbn3ZI{49 z@`tW#Xsn0w9IR4B*;a6oc8*JZ@m4tB`K-|vV^^5|c(M2Pvu=8rUu z?+$ADH64>@;q}Dv>vV8E3b@`KxSbogUD~>!x0f*Q3RgF<^r#$@GEuHo;)5%=MX(rpG~6 zSB$pq%XpyDt{?;R3bv#pM=w zy$1Zc)wrGlT<_O~n4OonUE2H{N@)>?6!HH457R1u|96Fiz{s z6pYgnExZ|L-|1nGrdycYP;?MvEt744#Mgpba;k=sUa&puPg`who!eCTRuMo+__G)qoVu3jZLBk6$4g83GFdCjKJ!_ z$=G!;3)2HvV|w8=%nn?Q*(FUzL;0egUljPYrey$Bi@wmJp`j#B`eE}B7RrCD7WkjD zTD0#Q{=r&GgR6ZVL};Mmq= znr`$ra9c$g2=$kfX_s7#18!+6@~m@HegdQoU$uFUFE7NN7bcQU4=e1N1c)saT2C2= z+>)7DRemjYIJw0m!4~!;bU3-y5mcw1P&S<0dNuz?&s^o<c-L8uL3)qPbJGrI zI0pu3{Gi*TwQH=$zwf$iKpME>=ZH|zdm2CJ_FS#{y85$62LefcL@FBGA1SZs7umW$ z<#eEInGY$iV$3R{M)ZSj^QyP=zSNVG+ev;*vi-4QL_g>@uNo90mP<6QCgl~n_<%Cu z$4^Lu7_aCV))G8u)xV>i+(TZT*h}+y6T%aEMz#Fr<^qrBKcDd9=Qr7izE8(^)*9{Yu|D(5NH2}eR^6>0*%7*(>94jph`d;tOL9b~d|T|u zj?nGk`}{E~cTX^!lv<4^`&JT;{5%6iO|6l^ED`)SLW%N^(9A_LfBF+MoZe#w0=hR| z0;}h~u{?A4I(UCHrMzKTEs%&`w!F#rCD=Bfx#E%BO8~?UW~(201LU{#dvA7m1NL3N z{&{cE8}M~T26ut&8!%;6B?z;pumh)=#U(Nv6)SF+5 zre*-ovWThK#hDLhS=sft^`|N4sh6^C$P1ee8%ZuPy>sN)v5-AjZDesGbRSMvdmeyC3 zP6R(xm-?$#o&ro^)_OPCl}KH5`++OFMhy1IP!|zoqWn(>JNuwZ?tJZ+AdfdvB{J(J zunp^e9n$m?T$>$nWSdI^@aucuHovM7*h$=P|L*q)Jls8jUFg*Vpx&<_q!CyFBo9A3 zm?rE2Hfk&RUAd75_RY(V$>z_8ORZSnT=(>Y!uoEiGAqhplW$pA(ZjnizPv}q)nY;1DA$aYMn{v?kg7$|6R$~YA}7}}KgeG+5uNk2aB}E- zd1Q@Ct=zmMT|{o$)vgVPC!t#xq{QUBmqWZ4-YvPFuojuyRrs!bD=%6m5jAV}QduPU zvK@D_tR8Y;w!ZfLlatUlr{?xuxcULGUT$LJJf??Oe!ncJ{&5Ogl&tpTjD{HUL3xeQ zjf%C%xN<|)+6jE?0$anWjoQ%1r5xydD%_BQC;gQd#ojn5%L?^0C}0eSXUd_pN8 z^`qEPiS^6S^l~#_rE>@}C3kXyW{(z^W@f8CBz@bzI)Fa70rsY z_F7&@oWhKJJ5$8pBNhmx9QB;4O_&iTRkMKI9xMp^<7p0Sw(W-@s)dp@Pd zxH*88pY+4f57zueODhF$JJYEj&)Ewd)TU`d`TvB;; zw>cZK!R$@$*q9^Guqbt}x(W+AF8S#$5ts1@&wF9@6CXpsb&d?Zch&3^L^pb^NHa6? z_7sx2OvVR(xgugWWyJuTo|j$sJfsU`224`W;y4REnXhp0mv%$e)J6xNwUwl$1z*o3 z>SrB>4!|ShYQ(@xlz*1PGJl-I%yt32w@Yl$11{R&n&ykPMe}V( zr-9FBXKoWhkGlC!6P-B#UZskiN-J|itxtU27`|{m>a_gD|)SwCME>8&(DO9l3+rIsy6 z^P{hvGwYa$wCbLDS>nAJbv{1!{pndM=&XQ5b)7I4q<_Jhii=;?qqWPX`x?GfMjvgR zr}aj;69``%*v{{;3Yq9#)yF5Tgf8!A?|HRoJTfJ65_?^=5t^%h=RULc3UqbIC9M@| zOsKC|;g^l!hRC#fWBE^kN~oLp&KR32X2exdGkB-nDpdQurL@USIrLcplS$?6G3Zte zd)aN%jF5}lA4fY={MF`q>cMj74OkNY>ih0&wzS=2SWwzy6hm>ipPUu3Soj-YxSt=k zAITpEL*Xv4MZ^$u<^giE7=aKx&aDF1Xo`J?q5DMcXv z`|7p2wZ>2m%_^8~ZwrN2B=96e+QYdTF(Qv!w!!3;YDf3WIKy-0-&=R)xxlI#V_0D5 z4C~FeR6lHXgR*Bnrz&lAg-SUZF&j;`!sbUIrA6FU@Y^#}-FkshN37OaWu%`a< z7g>j~aIAd0gW!A?@PF95^Kh!x_i_A$%wwdYGKXZ$5b9ZEDjAAIkuf22Man$SLoy^n z=6Mc@XKA7+^Q?%XP^d_R-|kbpPwlrl=eo}K{rk{$`e$8dANF?N>sj}D_FB(t-MbiR zB_p(s0)fLHW0^{h0f&%XIm32UJ1oV@jWrfvV{To5qU?fL%_0tIAjs2y@cQGKswc_CH*+d5Zxs8fNafCxbOW6#oC@&aHCEuKu5WjHiGPju~U zL^v#4G}_{@<0|AB@}3!;xdL-+9`~@gL_^LCMaz79!Xd)kP3UsO9p}4>-Q>pa*PDON z_R0P!yKFq>MA4apEhO# zyIWhLGd^Ylqm5Hu2URnGpZ1lcZ6#@7{6SysAaxQrc-lZuqA(tmO38IRlnMs0`xD1; z7YDG({p~DMu0GD#7I`1tkF6wMoXGb|vJ$bwrWm{ZzfQH0-C_zv%hL1%LXG;KIk(Bhj1i&BnNrX4PWotzuZ4wxjt z`TMo$H+H1L_q22uH1ZN*C? zk#XPvmtv^9auNu)f9$;7?kZ3%5iysgEdT=JzO~7i14rTcf?tkrKB(<3_MC{x0MoVH za*LiZI0IXT&~5+9*j$+ZA?GONRSEJ=Ka< znIxhE&JRKddq~XAIzbz42Az%lRT=IXKz=i&os^kJ07KR7&qAK6Kx(kLC>Z}D zxFw};p`+>!E?zpkljiU)oNVipT~b^*+iRm}7m|rPE!_?Pt>o7=Tg(mv*FYy-8vX;I zu1e7NsS5_@S@t;a7pMa~HQKFPBDBGa+K#!J{sW*brs$a{yCf(*YkfD#S_244I!@ZT z%EO*(WzC1e6ajMUqO02qyB_PkG)I1*|E>8hBM~AUbwZfz?VIRIr;} zc`14b^qnnlGO92GM%E(^r-d!Sp{5fW{RPH2mE^+d%35r;zbtgd7VQ;iAGUa4xBs+= z{L`*smnHbGP!QCDgRvjT`A!!A_JF(tSvhw%aM4dR5!tNc&cYo zX%(2F+0=F`>lS$B>qZ=>TMCM`lUyt@zh}Wr-^SCM2THk;FZO{f@amAyLbG5RxI)V5 z$dZ=~CJMZ_rmT+wM z;)}#Mm&ZTUm8`>t_T$x=UxPvH;ep-$Umqqn5PpS0eo`b^`~Kjx5j{ektWjSj$b|6T zJX7C%j0FikwfD{m5;kOrYu{5NVRl6IQtADILUv^Ax>1*KK06}wuGv~gjU9;)M4DeD zvLV)*1i2+tY)GDg>*5(oHst(i;y42(RwVLtnUJCZ3sSQ5^p1(dHe`F*dqaA{Ey$7m zlrc9xGay^j)B|E>sFC@%kFx5W2oc&|i7NG5GdMNrlS{_;TW3Gjq_x2bOT~>ur`Aj0 z459|DBx>v+`tfq~+Q5WTOFaDd%tdfwdVI}h2piLnJL(sd5PNMMyZxs}#Xr4IGyr&X z_SCT>Sx=)vbm1Vt(zkzK|Jwa+t>@wAlds!=uLpLpGW>AzUG-2=pU}04EExBE3vO*D z$AssIDd^3v8QV6_W$~n?H)>>d5EMQ)Hvbo z>QtgZVxQZ#A)8ncTY7@dkOe;&yfr#1vX&ajsLe}$Kgo);_6($S>6JrQjermaPAbHK z@o7!>5F29XC8{zbT?<3kx$@;mkt6FUMy>jb*%23GopaGqb@1@|BX2?yh!KU|!Bq+7 z9LNQg+o=^V8{ySM8?4zL5FiZww2EyXIgna?(^p3)nxPGMZ0s!UGMIGJ80**LM8Hm^ zDW|43NO(-G%ieYt-02VBwX=j1;ak$()6CZn*AeN{P465Bhdi#HQR(JH&Ww$4^yT(I zZbD0G&)F8RiN%oa;9E}QbD|N4Np=r>b>^n@x!D#-8{n|*%j+My_bVF&w86G=mHACO z$D#e~u<9Yq`<%0&cy?${D-^Z6yTD*M3+v1t*B!!qZl0B_iiA!htR44hu_RrF75=ro z(kD5QI1?o~W7m6-%Pr_xOAP@!xPD45c#Z?P!aTPjEwUQ6Rx<@(k0eH~kw4Z?J;Q-) z>Xl1pU@n5-{bRf!QF4?Z^L!^wDLY~-P3^#_9s%!Pc=+V+(`aeUD#YE*B_^=I_6tVqMIAQrXl4d6)l=p>Uo4QgPN9KXQ9iZfn+Y+Lu^{Bljx zN7aKAyRgOP+sW;#TVK6R>#D8qIOlRB=+9O2Ugw?Jl>=F|OV{;|aN%2J;XX+tOIbOPRXem;C7!s+NEqWB zBdh@pJSzvXYKN9jg3uGcL;G=t7BR~uVQe7ZX0l(a)@l}z;suUQyWZdmNTFGUgcU?pU3i~eK`J*^{_i=8ppPBmf{RTTuB&+tE56|tRe{Y(F6NzpP z?e}Mb;6$=&Px@8d+uqf91mGMdCf9d8$BwRV%cs7LuC-wyBh!`#3YLq2?#UVZ)~UN7 zO^`1nMyL~jdz{A)I1Pi7qIgfI1c!n0@Uo$&=P0oD-C}e9^e6~XSQ>2h9tD1j9_UV^ zQE=AtshH*75g_h=qU3DsOAtGkX}h)eB^bVBBzhqK1?a40tQ5aE3`{&rcza!jfYjE@ zVBO6=KrCes%ISK6K(mg)IL1N>D4ht};#>w|1;h5fR{0ML>FBj>|zg_7xKYv z+sr(__(JfY+0x_c;|#EudFl-EA`4JI)7|@uD;e0is#bU=#DUWK-0gv#J8;IfxJq~A z%Gq8U8uA^_cKx@=Fu-#)=x)xzINVJxE$4fNLAYmI_<`e3!Sk)eGvZYP;L=0|HMvbI z;JsBGDB0Q!?pIg8Ww5CSdrv8xKCW&LxDLhH65p->ds-feane`9EC(dT6tl?f>5Gkn zMe+|Im$(;eT1O+Kp4wVm(DxX6mnw)%$#lVeRnB-B7!7&%$&BX=MngUzcsl)fwhdl5 zsaG_E;%dkeo7(T!3cg**`C3WV28UD|-cV@PHw_)CyBD;pb|d_kS~mD;y#;AC?{wH) zAHyA+9`oZ}mO&&e?D!1XI1xG9w!YB1aoC^-E{*idAx4}Zd>i{DkPX!xh9e3~(Df*~ zFt9}iiC9;7dqiFlxqfqO(Zg#U%CtcikN^69#CT^R|D3cE6486mJW**q+N(Af&M_d4 zY;8A=_smm8tgd-}pl%^WbGF_sVH^@e_J{->m&#T{u66FBuB)a*xs-(r2?T@?;n`~$d#)WootR9^f5qEcTTvoEI;+!?|n-NU~ZeQjx5)0`!XwCp{0b((kt6i&z%LRATr zv4M`!FdIRv)mta8-I~L!vFelzz;H;skUX*|42LB2H`PkG`w0YRXLAqq38SgeuERmg zn~}{m`A!kFC`g*N+#UCh4gH`aTCU;oI}XYB^;QRVu>R{I|Mie-)I*x=dHMZb)Fs0qtxWm=j_Px0SgvNQae=H`2E34<-}W|M*#jA3a~dSkHU!Ev&v1eSa}o8I`c~{>{*?4?v^*hKfY#y$thnEDfYU({F|B)qS!uN! z<@WDC9o$WVa8A&eaz<&Q)sqb2RC^Irs7PT$d?6khCaDs`dszdCF``0JU8GU}H`y8o zZr~$K%}(Z@E-0cR<-L;l1iYvOXU33WEj~Ipoz=lcu8F7>eB4pZ&Wu{P1mi90d;&AA z${V%Aw2_jH+uDco`=Ob4{P}g`1PEJl9(^H`2IAXkaVS7o^&UMotLFAVPAeg(&UrX?xAR@ml*Z9wX));3yM?ebt0@kkICRjgz5e4^9jJ_RGNL$izFaC zPB9x6dIdD^O0Z5ZjRZST=60Rf5Wp|c*XLNI2wFajxe8ry1J#bUuTBuS!s!QZ<`M$; zL9{r5Wq(=#bPO?f-fW7&kaGi@H2lKh%V3g|6D(I@3Q>1$&f6>SO=pCfFa|}c;}1Q$ zv^N}*@*1~L9mYkGec|h?ev$6A0c;#SUO6$a90bSn%}JR=Ji(SXtXJ90q5y&JocmR> zRKRHAwqJK36&#qp^7*qxHt6qB*l=Ag2N*u@U4A!~1KO^tYPLq?0_x`Y7-*0K3b!l` z>jY(k<3{%=_qAq&yR}@j8y7M`6cJIz`Xd=Yq@*Np^ky1ham%=4K$irfYA(JQD~Shj z5TCRRqave39EkOe9RQsR<1*2}?^NWs^|9ZI$-gY}|22#B*VCq4d1$W@7y07OKHDp2 za$&+J>YArsH=$LQko>v23V6Tqwwba`9^5h8&S3}2paxaDwYXvuoE^>DTR@u%SuQ__ zY79%nAa}npwsm3PD(^d1+TGV6oh}De{nM-9jDG#c0irM{Zg;*l%Owunua{|gA(;fs zQ$1>UWsKru%t@A>`-bld0DuLYp07f(QB3 z#`6x)E$wcuD0%_D=po}tbGZaJ9C;1d8ZSbNC-sepgAZKpecUqK<_XI}WL~7I*+Keq z+T(f5#_-}6<;QR5*Dafare{(I~637}p9J~JB34Agq$cZ&` z1y6U&mBzjo!#Ph%9n(tu#pqreh5ND)JAHAP4-n8NFpcOI07&(G_cq(TuzTW2eC-SZ zWZBiP$e$I5r*d3O&P?n9Z)j*z^VJsjv7UG>e4gNSO~w7nIKlvRKtN@A?-P{*?7K!` zYoe+MV$(R8*WFNrEd(QPUWv;9p^J*n_n z$!p~_^2|V@1;@_Kh5;k4wUT5oPEf}{RD+}NQoA?* z2bA1=)bKk9{VCA3BPCTds+vae>cJH%`%MEFDaqP!RF4fKCGl?5zl3SyA|=HgFKWx>^}u}U(Fkvhl+-rXqxm8tA14g$630kMOH^m?S&Wn% znmMxl%j-gN3MCJN+MvOMBwl8Wlyr1q)c7JLi-;ZtrAW3y!k0wZn=w*y&&=7JFH-WE zLcnD*-A4FIUrfXfBPF+%`ec2PlFo^b?5o`GL1Kjuui7wDl8jE1_KTDh^t?2>2O}lV zj9*dvA|(kcIV`_O$qVY*TQ+Vkf|Z4|?*cGVQb4DaoE9S`)v4apyJMuJOF3W97b(d; zO=j^$N*2=9gc=1VfMLbN$9phRaw5hOp}|PWckRgK5{#7WEUI$&A|>~W+e$9}R!aUf z&s@_~vQVJ#`wix=UkUUb%+DIzoIk5?Z6l+(+FM4FZZRKv%7WW6a%fuB9(#HE*Fmk? zHSafAEVkfNYT!Vm@#DeKe;w4SUH*PKzp{N|FD?{&aUFc;^1s|utCv<1Ijrjl3BaY5 z+Z-ug|K*-iS+&EjebJ3`Slt$9xlni$zvW*CwQ7f-OQsutjo zs^4nnm)wOryhH_YkE*A5xe9+7)@r7ht@baJOSy3o`HvKH*Z)lXSJ+!?pJiqz-0%I% zy|IFzby5umur$c@5gT|e2H@-oKJpQ^RI(iwRDp{kgeHR zp#Uc+>m8JBt2#Ppw+UnKYNg>nDyIRLZ1r9(GTCoxZP7JXSI3ba*RmdZG`PboLdr zkNu6PR_$El?3HOL zeQtN16x!*U4YHjB2}R@Y0-lUh%Fo0*fivDj%^B}u(CK@9DrWaExG(E)?Hs1%Weh{| zUK#6AkWR4UoO8e^AQs+}R*Y$Znad#mt_s6G<8DRLAHcBB0zQ{{GE581{vFfEFs23O zP+8esL(gI0p24U6+zez0b8It`@sXOV?s`e8Z>oQrsi>4(WgQQ6p$ zQ3z7AiBpGR5ulK8TerTY5bT*rkX_tS0OUHJQzjPXfJ65c<*uZs0ff{!EtVh!h z-S~3A{6IYI0hk45w!R^_t(y&AkLgeao8|x)mty_vJ23q>&k&3~I+F&x?8{5I1u^Yb z&Reu7(Js4R`R0)k^DtBaE@PZ>sPQ$0gAHajh^_mACG(wY?VMN1?kKr2*nywzC3(5v> z-PTCb4yQeHMmU4JVA*XW^x?-gs2t|cuhfLgK|ghQ-NBCe9lxqJYlDWKm}%X^=|2Pa zG+kzjvK2&fDu&EHG`t5{Ph`rhgNNX$uHXF<5NAF|#z{scv@m$i^? zl|w!>hL;kY7DJ9?Ey|PhnTGTI zkH9zsdiPX{n(%fuBx%F$*!=5iNV&DP+(w!WXmfi|`@B6XN>2MR{(QA6vhQ?j3oE8m zW;73#x|G#0EFpWQ{^_kM(vz4#bBt|0Iw5v%s)n2&jXt(~fIajO@_0gTR`JymoF9Kz zGR!K1a$hnz>SA#a8P6l+AAj-z(q7cIU|SSMI|qqNJ~JpG&0FQK9xlp*m37ZbEj^{t zc)bLj7l&jJq7(02#_Th|sP4Et(Je7F!Ck7CLTVrK{>fu+mFhV#bKH;K{;D{tI&tm< z@|qWMz98p#se1vq(q>kqcnPBIT@i1S+UO9o;4rN_Ev>-afR`dGkr_Q6)f4R=g8OLd zA6=5fuLthmVxsO}TfX99xBuwbivPHAGz;0!(;QCP&o+N~h^`GfNUzT`P}F6Bx+v4n z9FyIFGVc0JQ+akfvOf$(YV>W%t?`YAL?2^n!vD)0HL{= zq82}J42>FK$gXMMi*9j!qZisY4f12hpP4c%A?WTC6b?jEXvy*VqDcY@BtDj*@PMf% zdQsbi`~tTOD#+x)-M1ef1&yCc*yA)1(_0<^_fAQn%;BD{+V1#BgM3s1<9h|Pk$aDH zk)s&8qr*Muk`D>WS$I+5e3d35pFzAqrHKbkauKEm~xjoiHJE?;h|pMNEm+61Ao$@! zc|2SOi1v=d#>)i)2;R42Ofdy?$b70?*Jrhs7Loa zJqnM9zczntWC|JbpiV)E86;Sgy~8qm1|C+vxS+^q4c|NT2xXPnz&B=*wd*h~7fE0H zzg$0P50wtwxNnbXx#+y^ZGQPN3phakjNO;o7}nU*@?A|h4uwsXpG%D&hRb3^_018a z@KPMhCOsiS@RadN#(>IE(9>Ype`Ed_kcbmG`(U>q&T&IAyBPR4QolCD%8gtPd;|4h zaK~PG$w4h1c(G1ZQyxhM`kt2hDZB~b3Zt+$7cu6)mGNasd#-@f#SdMk`XYhh_PRR# zh!7we^1LPxUk=PHSu_%exdDIVN|G6+E8Gum)MVa5AQ-mFC>RKUGqMu;bxz?h^oUgo z%kgj+u$PuZ1p!% zzc%!Wk9bj*B*kH{tFgJhW3w}8RcVjNIv5EGmrBMrNu+@J_2d0E=Ba?1!r$qEa5gY5 zTGv7>oCD5O$mwZ3%K;0|saPg23PqA&D?}1G0DEX&4aYr+L6O z8wV4b9(Sw1x(aUiR2+Y!8491oGmVs(UIQ)6(~0}rlfZLo)UVSa4&-ODX?zqd09E+q z6MplU%Onbf@!U4K;Em02`qHs9;7psQQhMYv&X6}6%5?o3sb3p~-u@ubn~qZfaA(!2 z3nw^EL2HrAvI4RWknvy{!sFr%KXDqA*xmJpr!Na$H@WT&2ldNKO{ct}e}R|YeOYf9 z)nKBPzTg3+EEHOd^L*i?15e=87Ejpr-pA_3s0$qIPF0Tfvx9N#@b_)3IR!^lwBm=N zj=?9!TDPT`^ud1ep{EMlPXeC29>8Ry^qI*CW1Q66^ov>ljnuD= zLjM;3nK22U2f#;85xtY%31auXXb5;E3NG*mJm)Gh0v!e2u3edjL2DWNI{I=Q5S{p< zS^maBz-naven3V7Sj32V;inw|`g;roqgxc<`ZL;j{;Z0C^W6E^+lSTR?#uIBZ~A$` z!tmq$XNE+8d1LGcGJhVxgNV8B@IL`6Hf9Pc=pF@p3P`9{hBer^p7s8%2s_|W@N{%4 zbOTQ6e{A5$@vB0=2FT*t?dHh)`>g-V7ypZV@t2pk|D}xoWy&b<`%(Yru*P3sj`*{P z<1eo%{mUKy``q!@b?i@}kH5Sf_NOq&U+%(xN#wsIausyp{Q1({k0kOd++X2t!XLmR ze?3k4!>Hshm)EUnHu?2&Fe?nrHBKnM;=2{%`V%F(r$)Fu_L+~pVnukl`EL8FgB+iPTvces=H2%*fDJ{(ZevyE#u0ZCv7p0 z!!XNRNtrjR2HK$4=!?L_ zFAXV=%60Hyxa0nJN`u=B+u{6ac_}fzanQk0@!$)yEH>Z!5YqQTSHgV=I%9~vv>%% zo?oB2nm24Jm&)k6qK`GVE4g0QT@)=~Vd@=C9sX zm5Oa5Ck!x8|NeOvS8Y4|z0&?Bs#UuZeevA=Hzvfm=#u^&SK&2S}L zANFdzYr|hhwR+iQ`{9tD6jNMw83}WW{zg=*b_51Z$R;1UVSsa7ZCrQeJ5l`gaiy|q zM_{k4nvn&sss5-zU^KHMXV^hiX@U&pnY>E-4m1N4F1 z@48^FbQSvR{#yI?aVyQDzbX@rTCtv`lUQv zeQDZa{(YInZ=zbY{L;U3ez;N{w{zyb5pdye?5kDlFAM5d+h|@-bjiR*wGy@e zT0vGZz`Sg0uPv5`3ou=GjLH3-uvYD?!z!n8l=^cT{CPH?mlMVMcfwk=(++RBInT=z~^7jzo}vJ6oRM=Scki0qT=a+`aBPxFb= zml%Wjl1Xg3X=fQQx*^^;=7ei7LvLGtpYi=B^=kuRUX$2mX-nn~-YV@qIFp+N4o+<> zlEL)h?64v@Q;+GxDetSC`1)okczV3n=YDz_*iDtnASGQ2LX|tFuW=OvWgnSaM~@VM znC)mFMO8W=5@Kl=_PYTT)B5H=Z7BpfP1YwoK^{n=9%#;8$_Kd@oZlR|UI={i8E+6a zWdLJJ>)Ce`S%A@REb%c!TAg=PW14cVKSM0tt^E`b^4xU_K?i^#ojDov*;YWU8JMKn-3<24 zQtFbs)&p2kmrCt_7MKp|_A5QB0Myd?yEzCdp}gwrIAu+5_J)7cCXa zH$ucFMKTJU+AcTb2C5@QZ|$EbaxtJMmftz<+Qx=h^P9>vC#fN~Lp1JWV>)ye!E+xk zo@GU!CA~Ttd0Q39K!c;R&k>{CdA9lwIEUe4?kRHd8C4`G;K1f>YWS#5a}-SzK0gYG z%?@q9d;Z>MFg$d!DrKEeGswhN_F%v>xAJ|bAkHP!ssLju{GSH zgp>+7y9>H^g6mT62|6xHqrUy*UJYt8$PTVH^@v>~K)*uRGNDKeU9MIw(*pYtx&xZ7 zUWIc&qde(^FNR_A8fl1fV;JVknF9=Y8|MLyW715TzaToCe0;&7j1GA++f46vs|_R> zid;wHnb852?FEtbzhjttR#;$#9l3vi!z4|~{eFq%*9Ym^fS9`PsbCr(tIyW4@%7J<5&$soImB0gaVC zwe!}=W!TBVl3|i$iY8FbKX?4dh1N^It*ft`0_QXIW;`ry(WAPrx#~LkP;I><2jfPb7ZyV(dV3Mai^q}!1Iw+=n|7L*@it?2{2sXS6gmMZx zVQ~>AK<=a#oFunGan9^NLd>si zbbr2&`}H9G5i6S7*c+ev{?!rdqOsf6<)eR$Zl%Q&{~6sHqnP7|%fqmFMd-&R5RXrC9A3^&=NawPhm}c#v~`#sNA%>z zUM%9Ka3D5-{dSBwG;*|dn|xso&4^CnIf~gsyjan5w=k_9XTs+{J`={YdQ`G$Ol{JG zk?kVlW0`7@elm(*F;oOveh4~u%s>QGpFF52w0S>x%szqi%N_-D9sPs{9kqa+R??)< zBCf^UGJo{@lSzJ;`?X<3ipr6IDr-2efTm%2QaE99{X zo7nDxzcjCmvE0_oXsSkz?Ev#b^lgIJ)foo_5W{9?&sG%ylsHT8SMV=(5B z+tHdu@vf%3Ru|vZGxn;LGsIIEy+Qt)g0n_E z<_5KgoO%;!aPzHP4xZ>T=v$O?IpKH(^riLFNgB$x#p7KXFuXe zzvHP;V6w=1*(MRncA+w?lVL!@Oc{CcJ`TnSSjZW_z6vxYMD1xTL*X?N-)R2B*MPh* z1DbI+2}EqNfTId=;H93b%2BQY@FFH;z-}xB%<)oJs$o24Z&TSmEt52W2MNfJDdT#~ zZ4Jo!zfo}3NUz8=$m5weKKo>Nlvfu3cm28&gmP~RwE-GNuO zkRxid`|*8-kVuRsq{W*D&?izlCoC9%DFbgOT{TOv{y93_Lu3y&C0`Z1>FxqzyB1ZS zB)WtCG)94wvnDumRHtIX%DqM`Dh6?*r-mROb8f$ic}AL`EUXM6m152j54oYM|Hr^u$&pMTqw% zO)Y6q0X%i<*bY}KI;;JpUCv)*r;Up_g?tg4&$ob1J@oRu&uH7z@D=onMdXD{?_J2+0 zf6-+A`fB&rJwXG9A9t|D~KHmTQ!*2i47h>Um{KPP_Cj-5CWE&Kb z^OgHFr(O#p>i16+6LxMy2=4=xt^!g-r@gCFu$3MW+QwE(QBQ_E-7fIya40c??Heh? z()>U85`Vz;W5?MK0K3J$j{Q08@9o(B^?g2ldtdD59RBvX*uSIZ+u!x`T|=5GYwy;y&Rj_`e*K!Y#(5^oFXz#Jg3nxfWZ{ne#_;Ru<=O$7 zgp@Rj7jHE~s(3i(8qp01$t|md{{v`N)Og7Mb<}G|Xi|^fM66xz0HM#XJA5!gbLS#~ z=^2dBT$jCv=H%-xY!(zd7`wk1mhKEVjh+Wis7bumJ7qH!CupTlEl zHT&{H5=LmIeLCIt03$T-pRpTr``8Q#m1E>LV}z!Oev!!+q3KB~tF-XE4KkYeRi&9%XX84rV8jL_U#!btn&eNsHaVjuIj z!O|@PWR)19Y4ZBig)c&rk9t`qeRC6ZNqp5wp z2+hfOs_HL7b6!hSSr;QTB~KnjzX;87asyHkjLsp?Uk^Ijt{3 z^TmNfM^iCE^LjG(%P&IH8-G^hi_qkaV|5WG%Y!#fhHfQdgyzn&<0p1vgl1g)Q@e1C z(4>yr;rK;p<_3co{=XHPzn;&oX`;!T^}RRkGUL|*>}2SuZ*O?+)N1dD`F2kb*X00t z97RU3ttSdwaaIm!)haumNpesmiTpf5~Z+D3rJr_b8H~>-p~7Rt#x1(@X6mTP{vr#qFC|DSk8oyO+M5 zZht-(ttEB&_3p}?+t2q6PB`~{zeT*_UR%xJQaj0~U*j$=n9bSFtd8AlE1qyyueIzp zjNT*fGJ?*sb2OxFMFu*Jh}|IuAR9zGTq zmVfk*Yyq$DbZIrN8?BGo@cpqrn^r34s~A^qP;oK8Rrmn@To(dKqW>zQRXc69iM(Dn zE|!fG(eSqFTi8MOi+!|eClQefCy$hJ+u=lXjABF@8_~Cie*b)WS|c4TEhQgw8gI=5 zah?Ht6&fITR2Eq z?j>O7+p45sF#^ba4e`CFUxE)UQXfAtjR2M1u8$3wMuDOndeEm7S{#bZnmTe=b3(= z^|b*h-Ou%e6l`;W3y-spzLdKG##|PuW6LwZ)|L*tAWX*+t!>KY#(Ft`<5ETP#Be^) zyMK`Fs;ukxGom>B{szz zeiMW<)~!61(tn@zwP7f2dDuUYF870QvAx0eF^|PFb+%1ms&x_PszW0Q=eD@+T)Mp@79* z6J|0`sCn2+hf1;@N{!lD@hUe%_gK#jqh2jgcF60BV@5xGcIk!_X?9Gb zhl;ir_kV&Vmp6vbiOC`gdKtZm50#K$9-{LPKYoT9;RnP7E94MyEw+RrW$&x9lJ2AeqE~XQ5`~r*!uGy(vYB8x)cLz2x0m1aol9srs@_JyK{4;bDBfSQ+HN zHirSN?K5zeDv<5kCULZvLq8+H6CiiimmhA=Z2-}9OZ|R&dr$@21xLZN0toYFIr(A? zwV1U#y4f%g`dj%S!~q9=EHy*u@O0sK65q3Ia-gdMcI zSDK6|=>IJ_d12kxxhHn}kEWi&LVPSD|7iY6W?+eD$39mAd6L4ZTSJME%>$;icc!Lb zVUFlV$7XNH5PEKF-Yx=U@w3b=hP!yky4@N)o`(X#uG#H#v|{TJ)rK&$JVpX!Ysot@ zpFt0Ja3JhrV&ycrbZ#A~xyE`Vi;qP zy*K_x9yUP6v-h7d`n&_uXDIle+OojsAA1Yy2%BJh%Ut=GSr1?k@VPsIQh{;#(dl^h zVn}6r8;?w*7|d;tb{b!7ukmd8NYrJBDfoG1SA(tw!yfcgxLz_DVDRidlOc}+0=?$K z;Pk7&eoCe-C$GGwW$=}@ZeKDur@vgd2h(2JyOl!D6SauDzK1eL>Va|o# z2@zKW!dbxkrrr8F2WtU|Sv^|Z*9y*^IoO}BDhDJ-RAIQ&9q@FK{``)n=YU;iX+y}E zB;@4URsSgXJ}7x!U}dyr8jvL<-DU0u;6oXgkidnAw=J91bB;_U#nUUZzFs>bS z`byt09F06SRVRawm>#2ehga_psV_h9+RKTDHm|oiAzMs_SVkQUYS|eAYA)q+C`b{Z zq!cI0h+^>(AMJusQ+5w{WhXkV=!S<9=_auDZJz}J*29->9nFRNW+T%T1V2Oet=rC1 zHdW%3uxw-_b9b<>G3xu$@npKwHuQn+v-m)S8Bq&mlF!rEZ+$#A#aS*@J_p&RH^{saL}s zZ#|xa4?6nzoh)_HvWWqTStoY1rzJWJN+3} z?SDP~JWL5mU#Hs`5x)y{YELwJx0C zcjAp7 z2WmzzlF+^L8EoC{QRm!Xf^cxg>6ZC!L0#XxJ`it6f{bThiZt7J0+Cr7y&W#Wgi34T z9oCnA0QRnbxAFE{6Qr@LG52jaDZ=JH%B!-04$*qbn(M-6h+H%wG^^Fwh-`~_dpNO| z8aXQ98YCI5f#}sNJ}8PNL}~I{yh{s6aV~fM$VtD+r=RB$e}a?!CX0R@*cv0UTCye7 z*h=$FF#H_bTG6p{x6kRQA36vr0+`p$@yLQ*nR_VU%Tn8G^tf~~O(7VxG_R;!SLC7g0ATO8eN4~6UQzZO!ngLcHk=PY6^ zp)t6d{K?Y{-h*U^)hkWl^UH}q@vF$;1U5J5o_!DHY!;#Qv z{bt49t|-{z)8!x&91f#6eKdDjdg7%0VUg|izexMqU{ZN^n6>fGo&X{X3F9}9dx0R6 z(h>jr(ZKnl%I8bkX&~SQu`~%crnBgrQ?2yt91zUg;OT6g3j#k+Uo*a)3l=QdBT6bI&LtI`zS7y? zpmH71vl^S(uw4W7lY*_{3_;*SLXrQLC+9&kopDeztv1eNee{#^d2HIhb)dj0GHZdu zR?h@^5gaIllT1g%^chG%!n zDYq0dl3&(C-z32FK`pbO`Wx_!bDoIKo&=boDLM+I!ol7rd<8SNuENF=LYhghSWx02 zwt(&og(-R*kL){Qz^>j$bDz@VfnBz##{;7%aFf7RfA~@kpkq*_bQ?(qCsk@hoU(I3 zVmG`n9ld(8or7%)QQS>Z~lw4uZ_Y!ADr9amf{X6RY|;`q??016n$xo9tMEy zjDU1jgf3jyMdE2$aRPF&5j~0-G=>j{KNrZeoq>_hE@`C}+Cm4lyxSj2&%$h#oP7>l zu5g%x$>O=CGaPUv(5M)9gihW2WC#+^!ci;xYBC~ws6}+`xFw?<%$wvec_VHD>jv)) zxtubEE${Eqzl}ZyIAd}sJ@L&!UhSq!8kbCQM!UYbruko_eQgx>TZ!WhVLh%OC@A>c z_79e@effk`NxBVOP#H??Y0!j{rB=N2XMmVc52WgxQUVgV3ENM-ZLeVEpZSU#n-F zP_E}>qG{9-;1&CMle01%s5^tdEi}&n)K{JkFKICZDm{`{G$|}`(oR$EE{IM0w}rGn z7xpzkUn?VX;x+Sb0;b4dMy;U~;G2Cp(5)&H1kue$nrEsfVs-K`32iETOHrB(&@x=wCuVy*(6%S~|Keq*N1)iQA0Cxk0pvlwtiPkgL=Rt$_oDk5HEIv&3*)t4(e77M1x zTW+dn;og|gD7jFA&G@%w-z#N&4XD?d1wq}IZ0en)6o!LXBKv9K7YJs~7!0Yl9 zrwyw>v-sRj;o>szDR&Q*eM~7h=wnDtN?HsKE-=OzOqy5jZt}9<& zSo_R%Z=u@v`|>MZ8(KTUby4c+9TLg)NYI76<4J~eNF|j9$Q<8{2&>Ruxk!v3aJ{Z0S> z{=CS#*Y^6JT!OVvP>Yv>7#g zf1Q1@S@fnBC!#!DP_ERuQ<9}awNNCySbe;pLXXSDy?|%6G z=Oy>7L)h^-kU$nbM~}|`ecf9@-Ciw({Xbq;8gh^Pw2A#c-sjyRJ8?#EHvSZE|q&H0yulA;`rXY@lHMkGu`>u|LYIW%-Xh6UQ}? zJdR=ckNv`s=(g9`g9TZSvmbxH?(k!`(Qln+em-UY0h)E?smLD&Syx;>x26s2mX+2C ztZ}}|v{Foef?TB!PbkBl{Z^bYt{t?hnbP2vQG{7LkpOL|pY;R5SAI8l4dHebT_R)s z4wP05cQY4w872ETfHOmiIGsZ1$ItT4rD5d#G#d^172y&>u+r#K%TYi%*dJLK(Z&N37TcP zVb}~iRSowt2p&DT2%9vJh{J@4*=h+IInAfj`D&4$uyC3q=os--zJ`PPs7mw;<-e;8hdTp}$+|fA~$@51O>jvbiiBs45>T2GRoENSQWXr&`CnCI0vVgg^I}L|do7jaTZz{3TH}%gy~bTJ*=h^>5$cNXwe@#>RuPFtvU0m1*1isrUVgGif7`um#n6CoU3qrIkWhYUb4!3V*KMg{rxqo-R&-WRDa(7 zSaOp{VAj4(Pb!pd)UY+_k6ljGRL<-g+IkM-?U1X-ch9B<9`hBnp9Amj_$AQs1M5WV zEKVgV=S?wGdQB2>Z35fnrDqZaS}q>*YvvR0TdKdw=$KC!9g{P-99l$(Yz+@hRCq#I zuYa{>_Sq*y;n2!}aL=cNQ;sOQsib!GxY8^W4V|aN?FRLqTyOAzd0a5doUHcIS@`>UL->)_c$s``15T#y@IYnr<^f8%{ z!YsVM7j0LMquc+E;?-_p`mq{X?nJYfrP)s^w-KMOh3xN@o=ns`p!wPNcq(zVXG+>+ zB9++Ye6i`gc?XC`XFs3qICnpx+}ZZQd;L`6an1YjmGS$CwNW8WTURB13U(!wm-y61S%q>kOW`kJu-9N@p3JN@%G2 z)Cw-xN~rEUJfW?}NwSHN7#G@?Lv#Fwc#EJ70H+et2Po!?zGSI*@gXq*V;r4Q$8-!ilqa@=y7l<5p z+ZT6rXA?D*14h*)j}mJ(`=zHys#k~4IM8pK%>wfA;Wp-aT`!S6ycTcn(?5fJJa2`Z zMrbA(eZji4-nLwFe^7!+)uMUwrn@5dy z=F_Jilg!^yzB;7CC*AMMipVGH9LbG~I*N2#9#wfAsUWhs+1(*!MJaJod^5mepp9r} zVQfXW^({q#HPn|aFD@huXL)=ZuHIK9`K9)NNx@d4Wi^~<7_@#s`pV_6&gozyikKPx zU`<9pk(P>y*P0Q}NiE~lv2Rm*i%fT!&VP5Vzv#lnYkdYa|3Dr!jkfq)-%`}QUe+j| zt`4HMDzB84UFF4d2WX8^KiWl<(n<6vaiN3AGF)Zqt%gcsD}xRtA4XpB;~6wgG@3;OUhSO8}&OhD@|FPtLY|Qzon2IkdRhx8#-8&)cD4p zpjx%XUZXzbY3a5XKO3;@-t-0oMZS6GYSwQdFaDg|CPLfQT&%gd%}eDM{YAD8mtFVm z`bd6HS1>r(!(7}rWaul4TehP5Ufm3irI(UR6|;<$GP{X;OMd#yextQ$z_Ru;MxDtg zGe$JZUsJz}cuD%q9?O>!BGrLMU+QaT6Cax@wP`rUQoQ>71Qmw>Q&Dl+#8MsSBI4u7 z#Qmz47UGD^Gvst@Em7&z?H$_OxI_{g4VxCb8i}t=(%Wq7RLOWPgI>PciLr%4joSa} zSO2SD{Y&dtqlR{qR%hdr^50Rw>b&T=$wMc1vC9r0odq!s#K}u?0vua67M)qJeDdqa zY2t3P+D+?Rs3)3uRV_Sag|4{d=H&0EC;E$xL(bhdR%$L9FeZMoSB0V2_4J_G8g8Cq zJ-Hg0Z40zSpH&p|w>s&FJLP!RaWwN4yY#u!`APQ%qJYY{t%(Dgi&I1GPo&Io7H`n6 znG$T=Ky*r}yF$_v4e=`z&8VnBW5hmM8w#BkloIoU)$58|X^1QLsb@S;au(a{*34e^ zT3tL~?!6DDI<>_I?N{d8ZRjUTyf!2`^1iNk@)ZB@k=<*G_P9)G+dxvldUn~mh50W! ziC0=oSv7sLlDK$=`^IgO0@k;)`_*ts=_EE+-B*94wvx!-Bg+2ys$*QLR8O>~f8CDXW)BtBY7=yOPs}&cb=T|s^72}Y z$Di~9)}KF>AboY^$E5sks9!kKZ&hx;&IKtvxulW&Sn}}x1HLul@+eT;%Ar2J?i-~c62%aG6qY!Jtm#9H{1&BAOH#yY`_BK>TuBk@ z*`lGY4k{w@(6N_y+t*)Ax(yt=)MaZd>2q!9_2@~ln^JU|o z^UQx`tpDF-tO*JAYI6+szod?J!pfb?&PFGZ+fwvWy5yvhVde9Wq&giX_h;QJd((3d zdCiUta9Dqc?5xsk{FSr~r1ANnR@#TQk{as%^+q<^K(?~F8m`lO8S%w;kII*0Ysgad zu1W^stBH~AgL*#_#gJCc2TPlitB8pi+R-Wz>xnNB;TzX=OdtYY*k!$+xSP1&aNjw3 zNhSV6zXe8hckU(*EKV)YKC*+*d1_Ga_9q#R)7P1a~@I?7$mgPb{ewBGlN&Scl7 zl^4b?=|x^h?K-eUi8-nKB0gq2(SgX2UpymM-IcI1Qn)XlIe};y8M^TEzOjs(>NQUn z=@$>O+*GIj9h19Bn~Af+>}C^h>sflXiJnI854>y;K72aq-6=-5i@g`w|Gas>jK`D7 z<%O?nT`e0)mM8b#pZ#n&>03G4Z@h;+xuzoa%d=hsiIX)-%FO#Zki7>y-8m=Mo*1|y z{^6|(JMwW{>yjEh`x0-TH&c$fJ&-sM?3kf9%9iM;<@CbL&4t+MUia9TkB)?=t)g1P zI^&7N6lYtr>thLGYPGMtE7_N<;)jef80-hboqa^tpCGW zR(z-8zhc(^aWU)9-(q4sUii`N|DI}AmPZZ$u5woVNzs4xtpEFZR+hQz-%`+uKRExl zRJ7tJ@4u4Pe@bKjD{1{vBK+r-wEp}ZOva7-Ur^J^vLgAH6}94T*8SC0t;`D0AL{g% zm$mAk3jS}XYqcTFKT;3Y@TKG5Q`mYytN-!a8*Y$em*0G-{4$p)n-{sQsZ3$(X+1+} zbr#;@zrV6|+VvwtjPkCLa~8TCQIk}*_ShBn`Iw}#wRPjzq>6c0$>MNPNUWr?b@?6L zO(P_gt%`NKTCEGZLHcBU)~_L{Y#q04xI!&Sg{Gs{aQ`tCS>(EvC+3Rx))wVhHgD%3 zspwotUes@ulub5@@XCLiQA0GoS?)lkfHC2_ZEt_?d9~@OcDmkmR|0t2m zzmttD8ynE} z_3?d4KVZIxGOFW3I(&utmVZS1xCcGZt= zRec)#{I6C2gvwX_@AMh|2zfhxmewj;pYR~2U@)?{WumAx?;pLeYpN-S+@476?ZQ3BiMRGZ19%DR|&^^PzksJ!4PWH8YCT`B`E}C*xzZ4_+=!3?-S7+j6UpHscjEs|yHIb02J{=04_IR>qBqNS79lCjVOa6=Z_?a`jeSNUeR`pU`s3gg%kD=s4{e{YFHeu56 z*W;M^(TTXvowJ$1sH1L#>TW}4XRU`*bn535UYUL$Q;p6pII zlb7gt{2tDmq$A6d-Qktt)oo@#@tT45p4OG&F<(nz!rT`Rv!`IN^VB44IASEVf}CT< zEg|j%&U4q7<#Bi3xJegSPVN!O)Wgh=yYt#Vd__?CNqS5A7b>sPg2`iS`!c(Z7~D9A zLP=t!f1&cDBv%Ve4c*RU0@-J-6eAqy<>@=sQ*sTvW0$YTj*mIiOy5eHUDd^=Ai4Ii z{IKmOW|CXXmdvN9@DeVPp^j(U?qNcc8|KDRT?uEJbV=qgReXDeNsJ$sHJ(Zn-lMeJ zFp1-2w+*Vp>4Uzkkip45rf1Tc4|@8FnxEkaZ zq`+VD1DKNFs*mJuBmE1NpCh>y!$J1zaTAlGx_D+JMS(+gj5(8nKh9wio8Db{j2aD) zXlKDB;*U$1#Hz1lqPh!?ug#s9L_}kmxWfZe<0u?FfkB;VxF1g+e|#kp6G68Ch_LG- zi$E8yX@}Zv$FtDtDvLlDu6>4Y^2HIwvItxk+ZHsJ&MUAXM=WKs7=$)N?&9GHJc|e` zSp?dU;!~ClCVhnWLW?dd2od#eYhi#BN(1DwwlA(OnQxgE~L3v6Y0!t_ycw~_vZ%HQsf-~{{Q_SP{=@pOj_l;O+P@iAw8 zy6#Cy{5^NqKEqJ|M%$`5mpp$c_K2uvP=0J%MA zdQ^%Qp5b>#1jsEyi8PDU2Oa6;Bn#xy$i13*d+_vj3`1z#9=Wu0RwNGm;*0>fJFv^x zbFXlq!w3Y(-GOxd(*qs2eWVOffJ(w8VEHH(7UR`e=A}n>NjQ-W7{fy3Et9Q5Yd5*D z@PyhVyqp)f%3yy>C*i!~;U+^B?vt=TPmGm83KU9sV^}p#h9}G@VS4`KWpDzOQu2^W z`WGsH!9#|}+bQ9|S??)B`aP+Hi;<8CGFTo~*+kmJs(XpOmkfiKR!W{XO8-LT4JOJ^ z*u8Sku#EWClHjTj9J7kvG8lehdF;D(9(6Z|$ak_x;bsZPncWnbEP1G9(fCm53R^?k zbycUJ`e*j=D5; zG#>te5ss>&jrl04U%{iFodu7MYGOg^%1G+y^`)n?>QY)70K;)b?4lHT8jl|t$cRVB zqeFVPd_3N0HX|M#kd<-wb*Tj!Tu+rT;<>fZgM6%vv#f()F5F6ZPwDY!oD~-gS=>6P z8obzonvvn|L&03Y;#R@k>KX5-okYpQS?OP>yjds%C0GNy?d`A`PjlA-z!s)}aGmoZ z9H4~jA6zJ3S;&g3R{d%lO|ha*FWgPr907=YwD0B3nICXiStJ9-Z7gT%Y@fZb-l_Ev3;qMLhl?Gal76hbOJ1k1XKvllCxI)Z05?*IeFZEH32>vy9jXE>YQJP2XJbj=jpdYT<(vIP5D(Ddu~ z@kCc2g>YfQ2A}X1ABQaA#s(J>Cr_~Rs#Uh(89tiSN!AB_8GjNY`8eAp---S>?a3)7 zjoU$<&S@yMBbZwkl1by15W2O&CFWu=gT|$8eh(d+;LW^uhKb_#Hrtf>YB;UUIVO#} z)7*#@%Uj892~anjF$P*}VaDLx?9m0ee984_i24o3TIJl*hodIFI z9S)9>otY5xdlCm%FOgZ0$;0BHgE`%T0A3mg_x@UBLj*sMg9Ee1ZOGs!ayt{^)u~%B z^=&TJWiX$wD+FuK(Ac~bZZ z58im6eq?l+X}0yDKOXG$fEJAF+nN#jP4Hl|TzW8?cF6W=()TN2hkXLUXxg{Br%K;< zX6)lUq#c>f+T=%Ump(&-;Y%OU!r7FKsQ*kF4sTi`@)+TKe9eE5pBdGSFdK_}8jelY zGDXWLQ{k{9ztO{SQyYj=fr252xYLFJUw~_(}?jzETS_9gz8vC-h@;dzq`sSm~2j zNUHgaPU4obrjh5h@tlSh(@ES$CZ3qQhni*K(r!x$jl`{EO)qS$P3@kw1KsrM1poDGuzY*bjii?Oa`Q z&q}xOAnMUafa30`t!bJ2@T}T@1}N^1ihP{dh(a|1S5oj78j3@;!n{gfQWUrpzwlLt z!r@u;S0CZT)##fHt6IsT7gn_ZZS(jpL*ipt^x`>q8CpwPFI|Z7)~k9=hIgd84bBm| z@-!Ylw@Nx_Hk`Wgz<%>lpuu>Em3diCG*F@eSxl<;;CazhG<*>K zh#t+VQSSyDTwg62)m0y;a>;dQM`w4avQjdc0>IH1sZ0a#6R1{+m!&n45H?MP1`}>i zaNL}!OM~%Hr_A@$hESahPtmG+KSfzembFnXLCh(=d}AD2gY;m2;P4x-VO;JdOn zBk|~R5RL28VPiw-?idORP6g0u7snRq@9wO?IyekJsRa%`6}^7gj~>|+ zhNHvK(@kqSb#sO5Ek|t-&h1j;2TWOvqqgV(6t_&xv*~#jN4?PnC~lhyS14JC-(87j zbd)fiY7Z7lOAR1Nuuo|wI%VQXDjL$!!h}j6Bnkw96YZ(rjW(mS#&3Yy8taa&{+c13NEbt%w;hG z+zQT#XDwu5yyXhMTcXiP7A4fK;3fGQ5PwUr;EZvtGd)G&eg*qZuPdDD@WA~3KzQ2x(E7oe4 zO9^;;{cUNnXwDiFDy{KYtA4atG-t8J{2o*+6at0*zhc?s>hanjmhFjav}iW9vT3tW zdTcoOR@(hKc2zVVvC8qUS&p}4+JIj)Hm!Oer>Rdx!*gdf@K-c0RlQpcV;QRogJ`j+ zQd!M#o`@g!3OyDTD)-%Ok5RF3S8MrT+HrB6ai?V~R*Cs*}>>%N!tU~Z8*+%%#)4jSc1%Ys{@{N^7D#6cH_(m>o2HLk`5zp)B( za#bI=QP6SNFA%pvSwwx?8Yp2@ecFQAFMl|lBV3>=G@9cHrj7aq7N$>dnLN(r7e}}} z!AVye3@&Co)p++Hj5%Z~p;7f1mGk$@p6(i!fmI{S{ z<6yJ*s}ju(pJPP{4SYba82Brk zU5J+04sS}u!$Fr%i|409`Z`OssdEfjmuCOUMYs&X@e&yPi^W5Q2HUB&r}`J3oI!K{ z0adr@bpFoS3Ou;cJbEy?%A`5)28*w@#)GZq(}U5pYwHJ@Q^9aF*H}OcW;3D=6@8^8fiV0AJ)Dmc!3lXy813kM zov7&B#{GC>W`xslY*Lh9k(fb+!}Cupq=n;F%+jSfDKjqMAl8z@|s zeFV`w+$bPS4<}8CqUXUaM_mTpi^oaX(R32G9Ub&NQ$fwHa7MI=p^=2?5qx3q^)dhw z>_;A!PmJ&s)8lDmVS)r7wEM55frJ|pyxsR&#R98VkxZv&mrV#NzW<|6yYdh z6F_lyRJY;In{d>=%>c#SQDHR6X##{ocPkAgK$YOKDRP?(MSv^8Deyz042!oc!5tva z9Wo@LwgfLjo1HY^Z|M@8BlLFBc*1=N_FHTc4JJ^S;56--Od|<1CYX}dZW>6SGQk)2 ziuck;yqyUShok#wjNg+cc&krHp#gbV6Vx8JR9ZAIZGu}-ZTHjS*}dt(QL=#YCOG<( z4$uJn#7S$|xTe&76T-^qFyZC|$IaF>T4p@dDb+N!j9SaVp&WR)isUxSd%C;t%;|1E zzSd6FZ7CI&4VG3PLC~NhY#?-f_+nJ_N<5?HN7+E=Qn6?wA$^YwX7qs#gsYZbfNBYz z(Y0eN8R75wExTzVEp>t!Z8^>cLYqd;3YAuC!R5)U6D%M$9eXIZr!Unwh#GW~1;t0r z;A>5Wr>a5m*}6*JkbWP}Qsp!YnoZeS*K{(bQ1B!avY~K&YuBgD2?re^S#rW_ubgI1 zl{PAWdx^%ac9^-xZ*uX$$pNLx}0Stp_g-Msi_7|s(X%=gbsSs zb$)TwqOKKKi9Kc|acf|=h*nBC)g_%RKW-%)pVng{PQ819jmoWupYDd8rS2idpa+Ir zVxe-YBK=$Ey8@{4RaU5AZ8XyO>WioTp6w{Y6cTO%?!Cd1x^SH|0e!OMCJ%76YPqUw zo)G#Y%~Y+ZT!LetWU^wpO=_g6`BWS`_7*Fa+p6gL32W^iWV2zpR8-@d z(;+<7p?6rB+*&%maMwH>d*d!Emb+tpb36Os*#0@JSniIcZ|uPr!p`4k!3t2T}Oum4HEnItHKhG;<0Si=M zcq&D@#~N zyi^*>tIu;57Q0rLUVbRu&V(a<_zM;oKgFINFjtXUV?gl0QWmgq-G<|`&TAGh5A}{c z{$TMU$rM%ffv5TfLalBK^VGD*RygU*8`=DRlF-$9gSW~#IB8WGD+yh#7r(fxiIb+j zWh3GG9Z&9-ZqdPh>h+G5gto4c!eps7Lz32eRubB}?_XnwQLTd)*K1Z1n}OGHjLXDX z7t7hQ;-ld34z;PGIu;)X&p(#)7|%2E0~?!7!d2z6TTv`{Cfq)&AWF5>fypZCU5N58Nph1d7(-JM{7qhi0Z zqtM|W6ZP>dzIdHb$%^83@tBiShT+(PZ>(5u8LxkeFvPKA

C{asHSd-MMXCeBC>j zy4|z_OU=9TY*=ABZV9MyitJRuK3-?ro+fzm!8O?qBTUHQwsU7CHmYzVhl{QTb@-60 zRq_`@j>J*}1S(EBE1Zv$H~lb>i+ za=Em;IOeexo@}8SJCxhY^{$F);apb@b}qM=&px1d4Cmgj&(7uUT>3^b{JKcz25ejb zY7Q5?*Be%k6yWCY8I@h*>cPC_+ziZZ>5^}kVs1A5n23LLflxb#Q$&gu8~wNR9L{X3 zo3g=$`#J3EMmlU{fr1Y2PM=t@!VDcI`#_g1S%HcUcb-@3vtfBVIvi=!4cLIcC+To~ z=hB>w&co87S$Ax~hT^5^`r!JiYsdy;_w;X9CWTVl-*Cu(V5RXB^`ocf_n{Uako?k! zEl=U54#(&2R%~P*s;-mN^gFe5hUF45t^dK!9v9%|KW>I+bvxuW^m3#Hj1c1127{mq zduvUY4W5ZvTMU9O-_LuwF2y18?GOmAv@fmFl>Ryed^@jKM+|~X z_7bI)ER=?wPzpMhAM2fx{#FOfbLVgQT3r6!L2Z(2MQEy&})kP zLvkmmFAW0xTF#oFYe7 zCT{_T!hU!tD^aNX!pnS}69VR~y**Oxk5eZJ&PgG|5Uz0Jh5h=`8G#CvU3mMrGaMmu zd+vsJ^NXp;9j1721OgPOxNyhFdo)7iZMSe3S&u+ ztCOinom#8IQ)=k}iPh~et66UDhf`EMQ3|?34w+pjeTx_7@(`z>D`ZbH@HU=H$^*EEgCuzrJyb9G$n|ja)B3&kvB$RGuRwCk1sgq6UyPE zu<&kjeG-;W6^D<*cDK9if#GR9?-*ow0A(vof1>kFq~9Nc2JybFt>vVoN0Q&9@4 zu~E~vJi{rieNYPOvHO~ADyJyW?!Vv^Tx1>O1Eibb5Oc#9O$Og7o}N_bhi9{88jeAS z_0{9MXH$y@xU`+&hhn&0_FlAIJskPMA4PJ@?8cgxZsW-DGf*VA%?^<}n@8Om;4Gdq z6GIBq*-k)bVlxms76y (body factory, solver kwargs) +CASES = { + "straight": (straight_body, {}), + "swept": (swept_body, {}), + "arc": (arc_body, {}), + "v3": (v3_body, {}), + "v3_viscous_noAT": ( + v3_body, + {"is_with_viscous_drag_correction": True, "is_with_attached_trailed_vortex_force": False}, + ), + "v3_yaw_rate": (lambda: v3_body(yaw_rate=0.4), {"reference_point": [0.1, 0.0, 0.5]}), +} + + +def numeric_outputs(results: dict) -> dict: + """Every numeric entry of a ``solve`` result, as float arrays.""" + out = {} + for key, value in results.items(): + try: + array = np.asarray(value, dtype=float) + except (TypeError, ValueError): + continue + out[key] = array + return out + + +def run_case(name: str, is_aoa_corrected: bool) -> dict: + factory, kwargs = CASES[name] + solver = Solver(allowed_error=1e-10, is_aoa_corrected=is_aoa_corrected, **kwargs) + results = solver.solve(factory()) + assert results["gamma_converged"] + return numeric_outputs(results) + + +def all_outputs() -> dict: + flat = {} + for name in CASES: + for corrected in (True, False): + for key, value in run_case(name, corrected).items(): + flat[f"{name}|{corrected}|{key}"] = value + return flat + + +if __name__ == "__main__": + np.savez(sys.argv[1], **all_outputs()) diff --git a/tests/BodyAerodynamics/test_panel_kernels.py b/tests/BodyAerodynamics/test_panel_kernels.py new file mode 100644 index 0000000..cc20943 --- /dev/null +++ b/tests/BodyAerodynamics/test_panel_kernels.py @@ -0,0 +1,303 @@ +"""The section force law as ``xp`` kernels (``VSM.core.panel_kernels``). + +Contract: + +* moving the force law into the kernels changed no output: every numeric + result of ``Solver.solve`` equals the one of VSM ``1201f6f`` (the commit + before the move) to round-off, for ``is_aoa_corrected`` True and False, on + a straight, a swept, a curved (arc) wing and the LEI V3 kite; +* the public quarter-chord induction matrices are the ones the force + directions are computed from; +* every kernel gives the same numbers with ``xp=casadi`` as with + ``xp=numpy``; +* a CasADi graph with the induction matrices as parameters reproduces the + per-panel forces of ``compute_results``. +""" + +import sys +from pathlib import Path + +import numpy as np +import pytest + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from panel_kernel_cases import CASES, REFERENCE_FILE, run_case, v3_body, arc_body # noqa: E402 +from VSM.core import panel_kernels # noqa: E402 +from VSM.core.Solver import Solver # noqa: E402 +from VSM.core.utils import thin_airfoil_strips, assemble_curved_bound_correction_AIC # noqa: E402 + +RTOL = 1e-12 + +_FORCE_KEYS = ("Fx", "Fy", "Fz") +_MOMENT_KEYS = ("Mx", "My", "Mz") + + +# ----------------------------------------------------------- no-change regression +@pytest.fixture(scope="module") +def reference(): + return np.load(REFERENCE_FILE) + + +@pytest.mark.parametrize("is_aoa_corrected", [True, False]) +@pytest.mark.parametrize("case", list(CASES)) +def test_solve_outputs_match_reference_commit(reference, case, is_aoa_corrected): + """Every numeric output equals VSM 1201f6f to round-off. Entries that + vanish by symmetry (side force of a symmetric wing, ...) are compared + against the case's force or moment scale instead of their own.""" + outputs = run_case(case, is_aoa_corrected) + prefix = f"{case}|{is_aoa_corrected}|" + expected_keys = {k[len(prefix):] for k in reference.files if k.startswith(prefix)} + assert set(outputs) == expected_keys + family_scale = { + "force": max(abs(float(reference[prefix + k])) for k in _FORCE_KEYS), + "moment": max(abs(float(reference[prefix + k])) for k in _MOMENT_KEYS), + "other": max(abs(float(reference[prefix + k])) for k in ("cl", "cd")), + } + for key, value in outputs.items(): + ref = reference[prefix + key] + assert value.shape == ref.shape, key + if not ref.size: + continue + if key[0] == "F" or key in ("lift", "drag", "side", "drag_induced_trefftz", "bridle_line_forces"): + family = "force" + elif key[0] == "M": + family = "moment" + else: + family = "other" + scale = max(np.max(np.abs(ref)), 1e-3 * family_scale[family]) + np.testing.assert_allclose(value, ref, rtol=0.0, atol=RTOL * scale, err_msg=key) + + +# ------------------------------------------------------ quarter-chord induction +@pytest.mark.parametrize("factory", [arc_body, v3_body]) +def test_quarter_chord_matrices_are_llt_plus_curved_correction(factory): + """The public matrices equal the trailed system on the quarter chord plus + the curved bound-vortex correction, assembled independently here.""" + body = factory() + va = np.array([p.va for p in body.panels]) + va_norm = np.linalg.norm(va, axis=1) + va_unit = va / va_norm[:, None] + core = 0.05 + public = body.compute_quarter_chord_AIC_matrices(core, va_norm, va_unit) + + llt = body.compute_AIC_matrices("LLT", core, va_norm, va_unit) + rows = lambda name: np.ascontiguousarray([getattr(p, name) for p in body.panels], dtype=float) # noqa: E731 + strip_x, strip_w = thin_airfoil_strips() + correction = assemble_curved_bound_correction_AIC( + rows("aerodynamic_center"), rows("bound_point_1"), rows("bound_point_2"), + rows("TE_point_1"), rows("TE_point_2"), rows("y_airf"), + np.array([p.chord for p in body.panels]), va_norm, + float(body._aerodynamic_center_location), core, strip_x, strip_w, + ) + for k in range(3): + assert public[k].shape == (body.n_panels, body.n_panels) + np.testing.assert_allclose(public[k], llt[k] + correction[k], rtol=1e-13, atol=1e-15) + if factory is arc_body: # curved: the correction is not negligible + assert np.max(np.abs(correction)) > 1e-3 * np.max(np.abs(llt)) + + +def test_quarter_chord_direction_angle_uses_the_public_matrices(): + body = v3_body() + solver = Solver(allowed_error=1e-10) + results = solver.solve(body) + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + gamma = results["gamma_distribution"] + v_qc = solver.va_array + np.column_stack([A @ gamma for A in AIC_qc]) + alpha_qc = np.arctan2( + np.sum(solver.x_airf_array * v_qc, axis=1), np.sum(solver.y_airf_array * v_qc, axis=1) + ) + np.testing.assert_allclose(np.ravel(results["alpha_at_ac"]), alpha_qc, rtol=0, atol=1e-14) + + +# --------------------------------------------------------- casadi == numpy +@pytest.fixture(scope="module") +def solved_v3(): + """A converged V3 solve and the arrays a force-law caller needs.""" + body = v3_body() + solver = Solver( + allowed_error=1e-10, + is_aoa_corrected=True, + is_with_attached_trailed_vortex_force=True, + ) + results = solver.solve(body) + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + return body, solver, results, AIC_qc + + +def _numpy_inputs(solver, results, AIC_qc): + gamma = np.asarray(results["gamma_distribution"], dtype=float) + aic = (solver.AIC_x, solver.AIC_y, solver.AIC_z) + v_rel = panel_kernels.quarter_chord_velocity(solver.va_array, gamma, *aic) + alpha, umag = panel_kernels.relative_flow( + v_rel, solver.x_airf_array, solver.y_airf_array, solver.z_airf_array + ) + panels = solver.panels + return { + "gamma": gamma, + "v_rel": v_rel, + "alpha": alpha, + "umag": umag, + "cl": np.array([p.compute_cl(a) for p, a in zip(panels, alpha)]), + "cd": np.array([p.compute_cd_cm(a)[0] for p, a in zip(panels, alpha)]), + "cm": np.array([p.compute_cd_cm(a)[1] for p, a in zip(panels, alpha)]), + "AIC_qc": AIC_qc, + } + + +def test_kernels_casadi_equals_numpy(solved_v3): + ca = pytest.importorskip("casadi") + body, solver, results, AIC_qc = solved_v3 + d = _numpy_inputs(solver, results, AIC_qc) + x, y, z = solver.x_airf_array, solver.y_airf_array, solver.z_airf_array + chord, width, rho = solver.chord_array, solver.width_array, solver.rho + + col = lambda v: ca.DM(np.asarray(v, dtype=float).reshape(-1, 1)) # noqa: E731 + mat = lambda a: ca.DM(np.asarray(a, dtype=float)) # noqa: E731 + + def same(ca_value, np_value): + ca_array = np.array(ca_value, dtype=float).reshape(np.shape(np_value)) + scale = max(np.max(np.abs(np_value)), 1e-300) + np.testing.assert_allclose(ca_array, np_value, rtol=0, atol=RTOL * scale) + + # flow + same(panel_kernels.relative_flow(mat(d["v_rel"]), mat(x), mat(y), mat(z), xp=ca)[0], d["alpha"]) + same(panel_kernels.relative_flow(mat(d["v_rel"]), mat(x), mat(y), mat(z), xp=ca)[1], d["umag"]) + v_qc_np = panel_kernels.quarter_chord_velocity(solver.va_array, d["gamma"], *AIC_qc) + v_qc_ca = panel_kernels.quarter_chord_velocity( + mat(solver.va_array), col(d["gamma"]), *[mat(A) for A in AIC_qc], xp=ca + ) + same(v_qc_ca, v_qc_np) + alpha_dir_np = panel_kernels.direction_angle(v_qc_np, x, y) + same(panel_kernels.direction_angle(v_qc_ca, mat(x), mat(y), xp=ca), alpha_dir_np) + + # force law, both direction laws + for alpha_dir in (None, alpha_dir_np): + np_out = panel_kernels.panel_forces( + d["alpha"], d["umag"], d["cl"], d["cd"], d["cm"], chord, width, x, y, z, rho, + alpha_dir=alpha_dir, + ) + ca_out = panel_kernels.panel_forces( + col(d["alpha"]), col(d["umag"]), col(d["cl"]), col(d["cd"]), col(d["cm"]), + col(chord), col(width), mat(x), mat(y), mat(z), rho, + alpha_dir=None if alpha_dir is None else col(alpha_dir), xp=ca, + ) + for ca_value, np_value in zip(ca_out, np_out): + same(ca_value, np_value) + np_vec = panel_kernels.section_force_vectors( + d["alpha"], d["umag"], d["cl"], d["cd"], d["cm"], chord, x, y, z, rho, + alpha_dir=alpha_dir, + ) + ca_vec = panel_kernels.section_force_vectors( + col(d["alpha"]), col(d["umag"]), col(d["cl"]), col(d["cd"]), col(d["cm"]), + col(chord), mat(x), mat(y), mat(z), rho, + alpha_dir=None if alpha_dir is None else col(alpha_dir), xp=ca, + ) + for ca_value, np_value in zip(ca_vec, np_vec): + same(ca_value, np_value) + + # attached trailed vortex force + legs = np.random.default_rng(0).normal(size=(body.n_panels + 1, 3)) + v_ind = np.random.default_rng(1).normal(size=(body.n_panels + 1, 3)) + same( + panel_kernels.attached_trailed_forces( + col(d["gamma"]), mat(solver.va_array), mat(v_ind), legs, rho, xp=ca + ), + panel_kernels.attached_trailed_forces(d["gamma"], solver.va_array, v_ind, legs, rho), + ) + + +# ------------------------------------------- casadi graph == compute_results +def _boundary_shares(n): + """(n, n+1) split of boundary loads onto panels, compute_results' rule.""" + shares = np.zeros((n, n + 1)) + shares[0, 0] = 1.0 + shares[n - 1, n] = 1.0 + for j in range(1, n): + shares[j - 1, j] = shares[j, j] = 0.5 + return shares + + +def _compute_results_at(body, solver, gamma): + """compute_results for this ``gamma`` with ``alpha``/``umag`` evaluated at + the same ``gamma``. (``solve`` passes the loop's last alpha, which lags the + returned gamma by one iteration, i.e. agrees with it only to the + convergence tolerance.)""" + alpha, umag, _ = solver.compute_aerodynamic_quantities(gamma) + va_norm = np.linalg.norm(solver.va_array, axis=1) + return body.compute_results( + gamma, solver.rho, solver.aerodynamic_model_type, solver.core_radius_fraction, + solver.mu, alpha, umag, solver.chord_array, solver.x_airf_array, + solver.y_airf_array, solver.z_airf_array, solver.va_array, va_norm, + solver.va_array / va_norm[:, None], solver.panels, False, + solver.is_with_viscous_drag_correction, solver.reference_point, + solver.is_aoa_corrected, + relative_velocity_array=solver.compute_relative_velocity(gamma), + is_with_attached_trailed_vortex_force=solver.is_with_attached_trailed_vortex_force, + ) + + +@pytest.mark.parametrize("is_aoa_corrected", [True, False]) +def test_casadi_graph_reproduces_compute_results_forces(is_aoa_corrected): + """Graph: gamma, va and the 3/4-chord and quarter-chord matrices are MX + parameters; ``v_rel`` and ``v_qc`` come from ``quarter_chord_velocity``, + the section forces from ``Solver.casadi_panel_force_function``. Evaluated + at a converged gamma it reproduces compute_results' panel forces and + moments (the attached-trailed part removed, as it is not in the + function).""" + ca = pytest.importorskip("casadi") + body = v3_body() + solver = Solver(allowed_error=1e-10, is_aoa_corrected=is_aoa_corrected) + gamma_converged = solver.solve(body)["gamma_distribution"] + results = _compute_results_at(body, solver, gamma_converged) + n = solver.n_panels + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + + gamma = ca.MX.sym("gamma", n) + va = ca.MX.sym("va", n, 3) + A_cp = [ca.MX.sym(f"A_cp_{k}", n, n) for k in range(3)] + A_qc = [ca.MX.sym(f"A_qc_{k}", n, n) for k in range(3)] + v_rel = panel_kernels.quarter_chord_velocity(va, gamma, *A_cp, xp=ca) + v_qc = panel_kernels.quarter_chord_velocity(va, gamma, *A_qc, xp=ca) if is_aoa_corrected else v_rel + forces_fn = solver.casadi_panel_force_function() + force, moment_local = forces_fn( + v_rel, v_qc, solver.x_airf_array, solver.y_airf_array, solver.z_airf_array, + solver.chord_array, solver.width_array, solver.rho, + )[:2] + graph = ca.Function("forces", [gamma, va, *A_cp, *A_qc], [force, moment_local]) + + force_num, moment_num = graph( + gamma_converged, solver.va_array, + solver.AIC_x, solver.AIC_y, solver.AIC_z, *AIC_qc, + ) + force_num = np.array(force_num) + moment_num = np.array(moment_num) + + shares = _boundary_shares(n) + f_at = np.array(results["F_attached_trailed_distribution"]) + assert f_at.shape == (n + 1, 3) # AT force on (the default) + expected_force = np.array(results["F_distribution"]) - shares @ f_at + reference_point = np.zeros(3) + ac = np.array([p.aerodynamic_center for p in solver.panels]) + at_points = body.compute_attached_trailed_vortex_forces( + gamma_converged, solver.rho, solver.core_radius_fraction, solver.va_array + )[1] + expected_moment = ( + np.array(results["M_distribution"]) + - shares @ np.cross(at_points - reference_point, f_at) + - np.cross(ac - reference_point, expected_force) + ) + f_scale = np.max(np.abs(expected_force)) + np.testing.assert_allclose(force_num, expected_force, rtol=0, atol=RTOL * f_scale) + np.testing.assert_allclose( + moment_num, expected_moment, rtol=0, atol=RTOL * np.max(np.abs(expected_moment)) + ) From a801b5ba31e2d53fa167118df04c008912d027dd Mon Sep 17 00:00:00 2001 From: ocayon Date: Tue, 3 Nov 2026 14:12:26 +0100 Subject: [PATCH 3/3] perf: vectorised Trefftz drag and projected area, cheaper Newton residual calls - compute_trefftz_plane_induced_drag: the O(n^2) per-panel Python loop with np.cross becomes one vectorised evaluation (1.3 -> 0.2 ms at 36 panels). It is a far-wake diagnostic computed on every solve; output unchanged. - Wing.compute_projected_area: vectorised over the strips (1.3 -> 0.1 ms). Also stops normalising the caller's z_plane_vector in place. - gamma_loop_casadi_newton: the inputs that are constant within a solve are converted to CasADi matrices once, and the loop calls a cached variant of the residual function with only the six outputs it reads (same graph). ~185 -> ~77 us per residual evaluation. Numerical, not physical. Sums are reordered, so outputs move at round-off only: every solve output still matches 1201f6f to <= 1.5e-14 relative. V3, 36 panels, casadi_newton: ~24 -> ~17 ms per solve. Generated with Claude Code --- src/VSM/core/BodyAerodynamics.py | 22 ++++++------- src/VSM/core/Solver.py | 53 ++++++++++++++++++++++++++------ src/VSM/core/WingGeometry.py | 28 ++++++++--------- 3 files changed, 66 insertions(+), 37 deletions(-) diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index 467cc95..a4d190e 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -1299,17 +1299,17 @@ def compute_trefftz_plane_induced_drag(self, gamma, rho): ).astype(float) d_gamma = np.concatenate(([gamma[0]], gamma[1:] - gamma[:-1], [-gamma[-1]])) q = origins - np.outer(origins @ e_inf, e_inf) # projected onto the plane - drag = 0.0 - for i in range(n): - mid = 0.5 * (q[i] + q[i + 1]) - d = mid - q - r2 = np.sum(d * d, axis=1) - v_t = np.sum( - (d_gamma / (2.0 * np.pi * r2))[:, None] * np.cross(e_inf, d), axis=0 - ) - # trace segment in the sense of the bound vortex (section i+1 -> i) - drag += gamma[i] * np.dot(np.cross(v_t, q[i] - q[i + 1]), e_inf) - return 0.5 * rho * drag + # velocity of all wake vortices k at every trace-segment midpoint i + mid = 0.5 * (q[:-1] + q[1:]) # (n, 3) + d = mid[:, None, :] - q[None, :, :] # (n, n+1, 3) + r2 = np.sum(d * d, axis=2) + v_t = np.sum( + (d_gamma[None, :] / (2.0 * np.pi * r2))[:, :, None] * np.cross(e_inf, d), + axis=1, + ) + # trace segment in the sense of the bound vortex (section i+1 -> i) + drag = np.sum(gamma[:n] * (np.cross(v_t, q[:-1] - q[1:]) @ e_inf)) + return 0.5 * rho * float(drag) def compute_attached_trailed_vortex_forces( self, gamma, rho, core_radius_fraction, va_array diff --git a/src/VSM/core/Solver.py b/src/VSM/core/Solver.py index 42d4342..96d3766 100644 --- a/src/VSM/core/Solver.py +++ b/src/VSM/core/Solver.py @@ -202,6 +202,7 @@ def __init__( # rejected and the step quartered (0.5 = residual may at most double). self._newton_reject_ratio = 0.5 self._casadi_newton_cache: dict = {} + self._casadi_newton_loop_cache: dict = {} #: Diagnostic: did the last casadi_newton solve fall back to base? self.last_fallback = False #: Diagnostic: residual evaluations used by the last casadi_newton solve. @@ -1359,6 +1360,31 @@ def casadi_panel_force_function(self): self._casadi_newton_cache[key] = fn return fn + def _casadi_newton_loop_function(self): + """The residual function restricted to the six outputs the Newton loop + reads (``h, dh_dvrel, mu, alpha, umag, gate``): the same graph, minus + the force-side lookups it would otherwise evaluate on every call. + Cached per polar set alongside the full function.""" + import casadi as ca + + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + key = self._casadi_newton_key(alpha_tables, cl_tables, cd_tables, cm_tables) + loop_fn = self._casadi_newton_loop_cache.get(key) + if loop_fn is None: + fn = self._casadi_newton_function() + inputs = fn.mx_in() + loop_fn = ca.Function( + "vsm_panel_residual_loop", + inputs, + fn.call(inputs)[:6], + fn.name_in(), + fn.name_out()[:6], + ) + if len(self._casadi_newton_loop_cache) >= 32: + self._casadi_newton_loop_cache.pop(next(iter(self._casadi_newton_loop_cache))) + self._casadi_newton_loop_cache[key] = loop_fn + return loop_fn + def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: """Newton solve of the circulation residual with an exact CasADi Jacobian, globalised by pseudo-transient continuation @@ -1407,7 +1433,9 @@ def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: fallback to decide. """ try: - fn = self._casadi_newton_function() + import casadi as ca + + fn = self._casadi_newton_loop_function() except ImportError as exc: # pragma: no cover - depends on optional dep raise ImportError( "gamma_loop_type='casadi_newton' requires the optional 'casadi' " @@ -1424,14 +1452,19 @@ def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: stall_angles = np.full(n, np.inf) stall_angles_neg = np.full(n, -np.inf) laplacian = np.zeros((n, n)) - params = ( - self.x_airf_array, - self.y_airf_array, - self.z_airf_array, - self.chord_array, - self.width_array, - stall_angles, - stall_angles_neg, + # Constant within the solve: converted to CasADi matrices once, which + # halves the cost of each residual evaluation. + params = tuple( + ca.DM(np.asarray(value, dtype=float)) + for value in ( + self.x_airf_array, + self.y_airf_array, + self.z_airf_array, + self.chord_array, + self.width_array, + stall_angles, + stall_angles_neg, + ) ) aic = (self.AIC_x, self.AIC_y, self.AIC_z) @@ -1440,7 +1473,7 @@ def evaluate(x): # v_rel = va + AIC gamma (see _build_casadi_newton_function). v_rel = self.va_array + np.column_stack([a @ x for a in aic]) l_gamma = laplacian @ x - h, dh_dvrel, mu, alpha, umag, gate = fn(v_rel, *params, l_gamma)[:6] + h, dh_dvrel, mu, alpha, umag, gate = fn(v_rel, *params, l_gamma) h = np.asarray(h, dtype=float).ravel() dh_dvrel = np.asarray(dh_dvrel, dtype=float) mu = np.asarray(mu, dtype=float).ravel() diff --git a/src/VSM/core/WingGeometry.py b/src/VSM/core/WingGeometry.py index ecc56c4..81fb1fd 100644 --- a/src/VSM/core/WingGeometry.py +++ b/src/VSM/core/WingGeometry.py @@ -627,22 +627,18 @@ def compute_projected_area(self, z_plane_vector=np.array([0.0, 0.0, 1.0])): """ sections = self.sections n = np.asarray(z_plane_vector, dtype=float) - n /= np.linalg.norm(n) - - def proj(P): - # orthogonal projection onto plane with normal n - return P - np.dot(P, n) * n - - S = 0.0 - for i in range(len(sections) - 1): - A = proj(sections[i].LE_point) - B = proj(sections[i].TE_point) - C = proj(sections[i + 1].TE_point) - D = proj(sections[i + 1].LE_point) - # two triangles: (A,B,C) and (A,C,D) - S += 0.5 * np.linalg.norm(np.cross(B - A, C - A)) - S += 0.5 * np.linalg.norm(np.cross(C - A, D - A)) - return S + n = n / np.linalg.norm(n) # not in place: never modify the caller's array + + LE = np.array([section.LE_point for section in sections], dtype=float) + TE = np.array([section.TE_point for section in sections], dtype=float) + # orthogonal projection onto the plane with normal n + LE = LE - np.outer(LE @ n, n) + TE = TE - np.outer(TE @ n, n) + A, B, C, D = LE[:-1], TE[:-1], TE[1:], LE[1:] + # two triangles per strip: (A,B,C) and (A,C,D) + area_ABC = 0.5 * np.linalg.norm(np.cross(B - A, C - A), axis=1) + area_ACD = 0.5 * np.linalg.norm(np.cross(C - A, D - A), axis=1) + return float(np.sum(np.column_stack([area_ABC, area_ACD]))) def compute_flat_area(self) -> float: """Calculate the total flat (undeformed) wing area in 3D space.