Walk on Decomposed Subdomains: A Hybrid Monte Carlo-Deterministic Solver for Elliptic PDEs

Elliptic partial differential equations are ubiquitous in graphics and engineering, but remain challenging to solve on complex or evolving geometries. Traditional discretization schemes (e.g., FEM/FDM) provide stable, globally coupled solutions but require heavy meshing or extreme refinement to accurately resolve geometric detail. In contrast, grid-free Monte Carlo methods (e.g., Walk on Spheres/Stars) adapt naturally to arbitrary geometry and offer massive parallelism, but rely on long random walks whose variance grows rapidly, particularly in the presence of Neumann boundaries, leading to slow convergence. We introduce a hybrid approach that combines the geometric flexibility of Monte Carlo estimation with deterministic global solves that do not introduce additional stochastic error. Our method decomposes the domain into simple, regular subdomains and uses Monte Carlo to estimate local first-passage solution operators (Poisson kernels), where walk lengths and variance are inherently controlled by the reduced spatial scale. These local operators are assembled into a sparse global system whose solution is obtained via a deterministic linear solve that exactly replaces simulating discrete random walks throughout the domain. This global solve trades stochastic variance for a fixed, resolution-dependent discretization bias, yielding stable and reusable solution operators. As a result, our method attains accurate, geometry-aware solutions even on coarse discretizations, and enables efficient solves and re-solves by computing and updating only the local operators affected by the geometry and its changes. We evaluate the approach on complex two-dimensional domains, benchmarking accuracy and convergence against standard grid-free and grid-based baselines, and demonstrate applications to microstructure simulation and flow-based path planning and streamline visualization.