WuWan: High-Performance Layered Elastic Half-Space Solver

WuWan is a high-performance C++/Python library for forward and inverse analysis of layered (up to 5-layer) elastic half-space pavement systems. Benchmarked against the Technical University of Denmark's ELLEA solver to a coefficient of determination of R² = 0.999999 across 5,000,000 evaluation points, it recovers layer moduli from a complete deflection basin in 5–50 ms, and provides a robust, uncertainty-aware D-optimal procedure for Falling Weight Deflectometer (FWD) sensor placement.

It is built for:

• Pavement Engineering — FWD data back-calculation and moduli estimation.
• Research — large-scale parametric studies, sensitivity analysis, and probabilistic inversion.
• Structural Assessment — real-time moduli estimation for quality control.

The solver couples a modern C++17 backend (Boost and Eigen) with a thin Python wrapper, using analytical gradients and structure-aware linear algebra to reach millisecond-level inversions.

WuWan moves beyond a conventional forward/inverse solver by treating measurement uncertainty as a first-class citizen: it quantifies how noise propagates into recovered moduli (Monte Carlo) and optimizes the sensor layout itself to maximize the information content of the deflection basin.

To establish numerical reliability, WuWan (v0.30) was benchmarked against the established ELLEA (v1.00) solver.

The validation study covered 500,000 distinct structural combinations, yielding a total of 5,000,000 evaluation points. Structural parameters were sampled within the following physical ranges to span a wide spectrum of pavement conditions:

The results demonstrate exceptional fidelity. As shown in Figure (a), the calculated deflections align almost perfectly with the line of equality, achieving a coefficient of determination of $R^2 = 0.999999$.

Figure (b) details the relative error distribution. The absolute difference between the solvers is negligible for the vast majority of cases:

The minor deviations (orange points) in the 0.5%–1.0% range are isolated to extreme stiffness contrasts — specifically, a very soft subgrade ($E_5 \approx 15$ MPa) paired with significantly stiffer upper layers. In these rare scenarios, numerical discrepancies tend to grow closer to the load application point. Even under these extreme conditions, the maximum error strictly remains below 1%; for all standard pavement structures, WuWan maintains a precision deviation of < 0.5%.

To verify that the inverse engine is mathematically consistent with the forward solver, the noise-free deflections from the same 500,000 randomly generated 5-layer structures were fed back into the back-calculation routine, and the recovered moduli were compared against the known "true" values.

Figure (a)–(e) plot predicted vs. true modulus for each layer. Every layer achieves a coefficient of determination of $R^2 = 1.0000$, with a Median Absolute Percentage Error (MdAPE) of:

Figure (f) shows the per-layer relative error distribution: Layers 1, 2, and 5 stay tightly centered on zero, while Layers 3 and 4 — which contribute the least to the surface deflection basin — show a wider (but still sub-0.001%) spread.

Figure (g) breaks down the residual magnitude across all 500,000 cases on a log scale:

In other words, 99.99% of the 500,000 structures are recovered to within floating-point precision. The handful of outlier cases (41 total) with larger residuals are confined to Layers 3 and 4, consistent with the known equifinality of layers that have weak influence on the measured surface deflection signal — the same effect observed in the noisy back-calculation example below.

Test Platform: MacBook Pro M4, single-threaded.

Note: The solver is optimized for large-scale batch processing in support of sensitivity analysis and probabilistic inversion.

The core algorithm solves the Layered Elastic Theory (LET) equations using the following numerical techniques:

1. Hankel Transform — the integral transform is reduced to algebraic equations using high-precision Gauss-Legendre quadrature.
2. System Solving — instead of a generic solver, WuWan employs a custom LU decomposition tailored to the sparse, banded structure of the 5-layer system matrices, reducing memory overhead and computation time.
3. Gradient Computation — the Jacobian matrix is computed by analytical derivation of the stiffness matrix, enabling precise sensitivity analysis without the computational overhead or truncation error of numerical differentiation.

The diagram below maps how a call from the GUI flows down into the C++ core and back, and how the uncertainty/optimization layer is built on top of the same core.

flowchart TD
    classDef entry fill:#0b2447,color:#fff,stroke:#0b2447,stroke-width:2px
    classDef core fill:#6a4c93,color:#fff,stroke:#6a4c93,stroke-width:2px
    classDef blue fill:#2f6fb0,color:#fff,stroke:#2f6fb0,stroke-width:2px
    classDef lightblue fill:#a9c9e6,color:#0b2447,stroke:#2f6fb0,stroke-width:1px
    classDef gray fill:#595959,color:#fff,stroke:#595959,stroke-width:1px
    classDef module fill:#3b4b59,color:#fff,stroke:#3b4b59,stroke-width:2px

    GUI["WuWanGUI.py<br/>Input: layered system (E, nu, h), load,<br/>evaluation points (r) with per-point ON/OFF mask,<br/>noise / search settings<br/>Output: deflections, recovered moduli,<br/>optimized sensor layout"]:::entry

    subgraph CORE["src/ - C++ Forward and Inverse Core"]
        direction TB
        TABLE["bessel_table.h<br/>Precomputed Bessel-zero lookup table"]:::gray
        PROC["processing_function.h + math_fun.h<br/>Gauss-Legendre quadrature, zero-segmented<br/>Bessel integration, coefficient / derivative kernels"]:::gray
        STRUCT["structures.h<br/>ModelParams, CalcBuffer, SimResults,<br/>BackCalcParams / Buffer / Result"]:::gray
        FWD["forward_main.cpp + interand_solver.cpp<br/>Calculation()<br/>Deflection and analytical Jacobian<br/>for the 5-layer half-space"]:::core
        FUNCTOR["inverse_functor.h + projected_lm.h<br/>Bound-projected Levenberg-Marquardt<br/>residual / Jacobian solver"]:::lightblue
        INV["inverse_main.cpp<br/>BackCalculation()<br/>Moduli recovery from a<br/>measured deflection basin"]:::blue

        TABLE --> PROC
        PROC --> STRUCT
        STRUCT --> FWD
        FWD --> FUNCTOR
        FUNCTOR --> INV
    end

    GUI -->|forward / back-calc call| CORE
    FWD -. output .-> GUI
    INV -. output .-> GUI

    subgraph OPTMOD["Uncertainty and Sensor-Optimization Module"]
        direction TB
        MC["montecarlo_main.cpp<br/>ParalleMonteCarlo()<br/>OpenMP-parallel resampling of recovered<br/>moduli under triangular noise"]:::module
        SLO["WuWan_pavement_slo.py<br/>optimize_sensor_layout()<br/>Differential Evolution over a sample-average<br/>approximation of the Fisher Information Matrix"]:::module
        COMPARE["run_monte_carlo_at()<br/>Re-runs Monte Carlo at the initial vs.<br/>optimized sensor layout for comparison"]:::module

        SLO --> COMPARE
    end

    INV --> MC
    FWD -->|Jacobian for FIM| SLO
    MC --> COMPARE
    GUI -->|MC / SLO call| OPTMOD
    OPTMOD -. output .-> PLOTS["Violin plots, DE convergence curve,<br/>layout and modulus comparison figures"]
    PLOTS -. rendered in .-> GUI

This example performs a forward calculation for a 5-layer pavement system using WuWanGUI, the desktop front-end built on top of the same WuWan C++ backend.

Launching WuWanGUI.py opens the main menu, from which the user picks one of three analysis modules: Forward Calculation, Back Calculation, or Sensor Location Optimization.

On the Forward Calculation page, the user fills in the editable (white) cells of the layered system table — modulus, Poisson's ratio, and thickness for each of the 5 layers, plus the applied stress, load radius, and sensor offsets ($r$) for each evaluation point — then clicks Compute!.

Each evaluation point also carries an ON/OFF toggle in the Use column. Points switched OFF (the row greys out) are skipped by the solver entirely — their deflection is reported as 0 — which is useful for probing how individual sensor positions contribute to the basin.

The deflection results (pink cells) are filled in instantly. Clicking Show Profile Plot renders the deflection basin inline, alongside the input table:

This example performs back-calculation (inversion) using the Back Calculation module of WuWanGUI, recovering layer moduli from a noise-free deflection basin.

On the Back Calculation page, the user enters the load (Stress = 0.95 MPa, Radius = 150 mm) and the measured deflection at each of the 10 evaluation points. Here the deflection basin is generated directly from a known set of "true" moduli ($E_{true} = [8000, 400, 300, 200, 100]$ MPa), with no measurement noise injected.

Every evaluation point has an ON/OFF toggle in the Use column. Points switched OFF — a malfunctioning geophone, an implausible reading — are excluded from the inversion exactly (their residual and Jacobian rows are removed from the least-squares system, not merely down-weighted), so back-calculation runs cleanly on any subset of 5 to 10 points. At least 5 points must stay ON, since 5 layer moduli are being identified; the GUI enforces this before dispatching the solver.

Next, the user provides the layer thicknesses, Poisson's ratios, and an initial modulus guess for the solver to start from — deliberately different from the true values, to test convergence. WuWanGUI also displays typical modulus ranges for common pavement layer types as a reference.

Clicking Run Single Calculation (no uncertainty / Monte Carlo sampling) recovers the best-fit elastic moduli from the deflection basin and renders the resulting layered profile.

Observation: Since the deflection basin contains no error, every layer — including the deeper, less-sensitive Layers 3 and 4 — is recovered essentially exactly from the deliberately off initial guess. This isolates the intrinsic accuracy of the WuWan solver itself, separate from the effects of measurement noise.

Building on the same Back Calculation module, this example quantifies how measurement uncertainty propagates into the recovered moduli. Instead of a single noise-free deflection basin, WuWanGUI repeatedly resamples randomized (triangular-distributed) noise on the layered system, the load, and the deflections, and re-runs the back-calculation for each trial.

In addition to the initial modulus guess, the user defines a modulus search range (lower/upper bound) per layer and a thickness noise (± mm) to be sampled for each Monte Carlo trial.

The user also sets the load (stress) noise level and, per sensor, a radial position noise (r ± mm) and deflection noise (± μm) — simulating realistic FWD measurement uncertainty.

Rows belonging to points toggled OFF in the Deflection Bowl table are mirrored here automatically: the row greys out, its noise entries are locked and set to 0, and the values are restored when the point is switched back ON — so the noise table always reflects exactly what enters the Monte Carlo sampling.

Clicking Run Monte Carlo repeats the back-calculation N = 1200 times, each with a freshly sampled noise realization, and plots the resulting distribution of recovered elastic moduli per layer as violin plots (95% CI, IQR, mean, and median).

[UR] Tail Risk: < 0.3 Excellent | 0.3–0.8 Acceptable | > 0.8 Poor. [RR] Core Spread: < 20% Excellent | 20–50% Acceptable | > 50% Poor.

Observation: The subgrade (Layer 5) is recovered with excellent core spread (RR = 3.11%) despite the injected noise, while Layer 3 — the layer with the least influence on the surface deflection basin — shows a poor core spread (RR = 75.97%) and the largest tail risk (UR = 4.39). This mirrors the equifinality seen in the half-million-case validation above: layers with weak sensitivity to surface deflections are inherently harder to pin down once measurement noise is introduced, even though the underlying solver itself is exact.

This example demonstrates the Sensor Location Optimization module, which searches for the FWD sensor radii ($r$) that maximize the information content of the deflection basin for back-calculation, under realistic uncertainty in the layered system, load, and measurement noise.

As with the Back Calculation module, the user starts from a measured deflection basin (Stress = 0.95 MPa, Radius = 150 mm) and the corresponding layered system (initial moduli, Poisson's ratios, thicknesses). An optional "True Modulus" row lets the user compare the optimization against known reference values.

The user sets how many sensors are fixed (kept at their original positions) vs. free to move, the allowable search range $[r_{min}, r_{max}]$, the minimum spacing between sensors, and the Differential Evolution (DE) solver settings (population size, iterations, tolerance, random seed).

Sensor masking carries through here as well: sensors toggled OFF take no part in the optimization — they contribute nothing to the Fisher Information Matrix, occupy no dimension of the DE search, and reserve no minimum-gap spacing — and Number of Fixed Sensors counts enabled sensors only. Disabled sensors are reported with a greyed-out Disabled status in the result table and layout plots, keeping their entered position untouched.

As in the Monte Carlo back-calculation, the modulus search range, thickness noise, load noise, and per-sensor deflection noise are defined as triangular-distributed priors. These define the uncertainty that the optimizer is made robust against.

Clicking Run Sensor Location Optimization performs a robust D-optimal search (Differential Evolution over a Sample Average Approximation of the expected Fisher Information) that keeps the first 3 sensors fixed and repositions the remaining 7 within $[300, 3000]$ mm, subject to a 100 mm minimum gap.

This corresponds to a D-efficiency of 1.86× relative to the initial layout, shrinking the 95% confidence-ellipsoid volume of the recovered moduli to 21.2% of its original size.

Observation: The optimizer pushes the free sensors outward, toward the edge of the allowed search range, because the deepest and least-sensitive layers (Layers 3 and 4 — see the validation and Monte Carlo sections above) are best identified by sensors farther from the load, where their relative contribution to the deflection basin is largest.

WuWanGUI's Preview tab provides three additional diagnostic plots once the optimization finishes:

DE Convergence Curve — tracks the D-efficiency of the candidate layout (relative to the initial layout) over the 200 DE iterations.

Modulus Distribution Comparison — re-runs the Monte Carlo back-calculation under the initial and optimized sensor layouts side by side, showing how the optimized layout narrows the recovered modulus distribution for every layer.

Sensor Layout Comparison — visualizes the initial (top) vs. optimized (bottom) sensor positions along the radial axis. Sensors toggled OFF appear as greyed-out ✕ markers on the initial row only, since they take no part in the optimized layout.

• High-Performance Core — originally a Cython project, now fully rewritten in C++ for maximum efficiency.
• Analytical Gradients — analytical derivatives for Jacobian calculations, outperforming finite-difference methods in stability and speed.
• Advanced Back-Calculation — Eigen combined with high-speed C++ gradient providers, solving inverse problems in tens of milliseconds.
• Robust Error Modeling — noise injection for thickness, deflection, load, and sensor positioning to simulate real-world measurement uncertainty.
• Selective Point Masking — per-point ON/OFF toggles exclude individual sensors (a faulty geophone, an outlier reading) from the analysis exactly, consistently across forward calculation, back-calculation, Monte Carlo, and sensor optimization; any subset of 5–10 points is supported.
• Sensor Location Optimization — robust, uncertainty-aware Differential Evolution search for the FWD sensor layout that maximizes the information content of the deflection basin.

• Operating Systems: Linux, macOS, Windows
• Python: 3.8 or higher
• C++ Compiler: supporting C++17 (GCC 7+, Clang 5+, MSVC 2017+)

• C++ Compiler supporting C++17
• Boost Math Library
• Eigen3 Linear Algebra Library
• Python 3.x

git clone https://github.com/lewiswan/WuWan.git
cd WuWan

This method automatically sets up a build environment, downloads the necessary C++ libraries (Eigen & Boost), and compiles the project.

pip install .

Note: The first installation may take a few minutes as it downloads the Boost C++ headers.

If you are modifying the C++ code or reinstalling frequently, use this method. It disables build isolation to persist the CMake cache, preventing a re-download of Boost/Eigen on every build and reducing compile time to seconds.

1. Install build tools (one-time setup):

pip install cmake ninja pybind11

2. Fast install command:

pip install . --no-build-isolation --no-deps --force-reinstall

For comprehensive guidance on using WuWan:

• API Reference: detailed documentation of all classes and methods.
• Tutorials: step-by-step examples in the examples/ directory.
• Theory: mathematical derivations and implementation details in docs/theory.pdf.
• FAQ: common questions and troubleshooting tips.

• C++ Core Rewrite — transformed from Cython to C++ with Eigen/Boost.
• Forward Calculation & Analytical Gradients — high-speed forward modeling and derivative calculation.
• Deterministic Back-calculation — fast inverse analysis for moduli estimation.
• Monte Carlo Uncertainty Analysis — triangular-distributed noise injection on thickness, load, sensor position, and deflection, with UR/RR risk-spread reporting per layer.
• Sensor Location Optimization — robust D-optimal sensor placement via Differential Evolution over a Sample Average Approximation of the Fisher Information Matrix. This is a preliminary release of the module — the search heuristics and robustness criteria are still being refined.
• Selective Point Masking — per-point ON/OFF exclusion of evaluation points, applied exactly (residual/Jacobian rows removed, FIM/DE dimensions reduced) and consistently across all four analysis modules, with a ≥ 5-point guard for the 5 unknown moduli.

This project is licensed under the Apache License 2.0 — see the LICENSE file for details.

You are free to use, modify, and distribute this software for both commercial and non-commercial purposes, subject to the terms of the license.

WuWan builds upon decades of research in layered elastic theory and pavement mechanics:

• Libraries: built with Eigen, Boost, and pybind11.
• Validation: benchmarked against ELLEA by the Technical University of Denmark.
• Theory: inspired by foundational work from:
  • Levenberg, E. (2020) — Pavement Mechanics: Lecture Notes
  • Huang, Y.H. (2004) — Pavement Analysis and Design
  • Ullidtz, P. (1998) — Modelling Flexible Pavement Response and Performance

Special thanks to all contributors and the pavement engineering research community.

If you find WuWan useful, please consider giving it a star on GitHub — it helps others discover the project and motivates continued development.