Source-linked AI summary
CraftSPH: A high-accuracy and composable differentiable SPH solver implemented in PyTorch
Gen Matono, Shujiro Fujioka, Mayuko Nishio
TL;DR
Differentiable SPH solvers have advanced, but unified frameworks remain limited in their support for composable high-accuracy and implicit schemes. CraftSPH combines independent PyTorch modules into differentiable solvers, with forward benchmarks and inverse analyses demonstrating its applicability.
Problem
Existing differentiable SPH solvers offer limited flexibility for incorporating high-accuracy discretizations and implicit schemes, while their tightly integrated operations hinder customized solver construction.
Source: §1 · p. 3Method
CraftSPH implements SPH operations as independent PyTorch modules that users combine into solvers, with automatic differentiation supported across modules.
Source: §1 · pp. 3–4Results
CraftSPH reproduced benchmark behaviors across internal, free-surface, multiphase, and three-dimensional flows; inverse estimates had viscosity errors of 0.19% and 12.13% for WCSPH and ISPH, and Reynolds-number errors of 0.47% and 2.03%, respectively.
Source: §6 · pp. 36–37Takeaways & Limitations
Different SPH solvers can be constructed and reused without fixing the overall numerical procedure in advance, including extending the two-dimensional dam-break solver to three dimensions without major procedure changes.
Source: §6 · p. 36Validation mainly focuses on incompressible flows; broader higher-order operators, tensor-valued formulations, and coupling with other numerical methods remain to be developed and validated.
Source: §6 · p. 37
Abstract
from arXiv · showhide
Smoothed Particle Hydrodynamics (SPH) is well suited to a range of problems, particularly those involving large deformations in fluid dynamics. In recent years, in addition to the advancement of SPH formulations, the development of differentiable solvers has also progressed. However, unified frameworks that flexibly accommodate diverse numerical schemes, including advanced and implicit methods, while supporting continuous extension and updating remain limited. In this study, a high-accuracy and composable differentiable SPH solver, CraftSPH, has been developed. In CraftSPH, major computational operations are implemented as independent modules, which can be combined according to the intended purpose of constructing a solver. This design enables different computational schemes to be intuitively constructed from combinations of common components and allows newly developed high-accuracy discretization methods to be flexibly incorporated and extended. Furthermore, each module supports automatic differentiation, enabling applications to physical-parameter estimation and to optimization that combines SPH solvers with deep learning models. Accordingly, CraftSPH provides an SPH computational framework in which high-accuracy discretization, explicit and implicit computations, and automatic differentiation can be handled in a unified and extensible manner. To demonstrate the applicability of CraftSPH to a wide range of problems, forward analyses are conducted for Poiseuille flow, two- and three-dimensional dam-break problems, and a rising-bubble problem. In addition, the performance of automatic differentiation is evaluated through inverse problems involving parameter estimation for the Taylor-Green vortex and lid-driven cavity flow.
1. INTRODUCTION
CraftSPH addresses the limited flexibility of existing differentiable SPH solvers with independently composable modules for constructing advanced schemes. Its PyTorch implementation supports automatic differentiation, deep-learning integration, and forward and inverse analyses.
1. INTRODUCTION
- Parameter estimation and differentiable deep-learning models require gradients of numerical results with respect to model parameters.Differentiable programming represents solution steps in a computational graph, allowing automatic differentiation to compute those gradients.Source: §1 · p. 2
- Existing differentiable SPH solvers mainly use standard discretizations and explicit integration, with limited support for high-accuracy and implicit schemes and tightly integrated operations.These designs make it difficult to rearrange components into customized schemes and couple solvers with deep-learning models.Source: §1 · p. 3
- CraftSPH implements key SPH operations as independent modules that users combine to construct solvers and incorporate advanced numerical schemes.Particle variables and material properties are managed through a unified particle object.Source: §1 · pp. 3–4
- Implementing SPH and deep-learning models in PyTorch supports hybrid models, gradient computation, and comparisons under consistent execution conditions.CraftSPH supports CPU and GPU execution.Source: §1 · p. 4
- The paper evaluates forward analyses by combining modules into advanced solvers and inverse analyses by assessing automatic-differentiation gradients.The stated focus is fluid dynamics.Source: §1 · p. 4
- CraftSPH supports high-accuracy discretization, stabilization techniques, and both explicit and implicit schemes.The contributions also include a PyTorch-based differentiable solver for inverse analysis and deep-learning integration.Source: §1 · p. 4
2. DETAIL OF SPH METHOD
The section presents governing equations and spatial and temporal discretization for incompressible-flow validation. CraftSPH is a general-purpose framework whose components can also be combined for other governing equations.
2. DETAIL OF SPH METHOD
- The section describes governing equations and spatial and temporal discretization for the incompressible fluid flows used in validation.Source: §2 · p. 4
- CraftSPH is a general-purpose differentiable SPH framework in which users assemble components into solvers or integrate SPH computations into deep-learning models.Source: §2 · pp. 4–5
- Solvers for other governing equations can be constructed by selecting and combining components in different ways.Source: §2 · p. 5
2.1. Governing equation for fluid dynamics
For incompressible-flow SPH analysis, the governing equations are the Lagrangian Navier–Stokes equation and the continuity equation. The passage identifies the physical variables and external-force contributions used in those equations.
2.1. Governing equation for fluid dynamics
- The governing equations for incompressible-flow SPH analysis are the Lagrangian Navier–Stokes equation and the continuity equation.Source: §2.1 · p. 5
- The passage identifies time, velocity, pressure, density, and kinetic viscosity as variables in the governing equations.Source: §2.1 · p. 5
- External force includes effects such as gravity force and surface-tension force.Source: §2.1 · p. 5
2.2. Spatial discretization in SPH method
The section contrasts classical SPH, CSPH, and LSSPH as increasingly corrective approaches to spatial discretization, with accuracy and stability depending on particle arrangement and neighbor availability.
2.2.1. Classical SPH
- Classical SPH is simple, but irregular particle distributions reduce gradient convergence to zeroth order and Laplacian convergence below zero, potentially causing critical numerical errors.Particle motion can produce such nonuniform distributions during simulations.Source: §2.2.1 · pp. 6–7
2.2.2. Corrected SPH (CSPH)
- CSPH corrects the kernel gradient using a moment matrix, achieving first-order gradient convergence for disordered particles but not improving Laplacian convergence or exceeding first-order accuracy.The moment matrix can become singular near free surfaces or strongly biased particle distributions.Source: §2.2.2 · p. 7
2.2.3. Least squares SPH (LSSPH)
- LSSPH uses Taylor reconstruction and moving least squares to approximate derivatives of arbitrary order and accuracy in principle.Its derivative coefficients are determined through a weighted least-squares formulation.Source: §2.2.3 · p. 8
- In practice, LSSPH accuracy and stability are constrained by neighbor counts and numerical stability, so lower-order models and stabilization or model switching may be needed for distorted distributions.LSSPH requires more neighbors than CSPH and is more susceptible to instability; Type-A is more common for its smaller moment matrix and lower cost.Source: §2.2.3 · pp. 9–10
2.3. Temporal discretization of SPH method
The section compares explicit WCSPH and semi-implicit ISPH for incompressible-flow analysis. Their time-step and per-step costs differ because WCSPH is sound-speed constrained, while ISPH solves a pressure Poisson equation.
2.3. Temporal discretization of SPH method
- Explicit methods are cheaper per step but often stability-limited, whereas implicit methods require iterative solves but can permit larger time steps and lower total cost depending on the problem.The framework can also combine modules to construct methods beyond WCSPH and ISPH, including higher-order time integration and fully implicit schemes.Source: §2.3 · p. 10
2.3.1. WCSPH
- WCSPH closes its equations with an equation of state and advances explicitly, but its CFL time step is severely restricted by the artificial speed of sound.The artificial sound speed is introduced to limit density variations to approximately 1%.Source: §2.3.1 · pp. 11–12
2.3.2. ISPH
- ISPH uses a predictor, pressure Poisson solve, and corrector; its incompressibility-based pressure removes the acoustic time-step restriction but requires a global solve each step.The larger allowable steps may offset added linear-solver cost, depending on the problem, resolution, and solver performance.Source: §2.3.2 · p. 13
3. IMPLEMENTATION DETAILS
CraftSPH is a PyTorch SPH library built from independent modules rather than a monolithic solver with fixed equations or time integration. Users combine modules to construct solvers for their analysis.
3. IMPLEMENTATION DETAILS
- CraftSPH's core design is a library of independent components that users combine instead of relying on a fixed, monolithic solver.The components can be selected according to the analysis purpose.Source: §3 · pp. 13–14
- The modular structure is intended to support constructing different SPH solvers and incorporating new numerical methods.Individual numerical-method formulations are described in references cited in the corresponding subsections.Source: §3 · pp. 13–14
- Its modules cover neighbor search, kernels, differential operators, pressure computation, stabilization, boundary treatments, time integration, and output.The library also manages variables and particle types.Source: §3 · p. 14
3.1. Management of Information and Program Execution
CraftSPH separates analysis settings, particle information, and solver execution through configuration files and a particle object. This organization supports changing problem settings and numerical methods with relatively minor modifications.
3.1. Management of Information and Program Execution
- Configuration files specify analysis conditions and initialize an SPH_Particle object passed to a user-defined solver, separating settings from the time-integration algorithm.Required configuration items and consistency across particle data, solvers, and output settings are checked at initialization.Source: §3.1 · p. 14, §3.1.3 · p. 17
3.1.1. Variable management
- The particle object registers named physical quantities and stores simulation state, including particle count, time, step size, and neighbor-search information.Scalar and vector quantities can be defined with their tensor ranks and data types; reserved names prevent conflicts with internal variables.Source: §3.1.1 · pp. 14–15
3.1.2. Particle type (object) management
- Particle types and sets are configured separately from variable names, with automatic partitioning and grouping that avoids manual tensor concatenation or splitting.Users specify target particle types while configuration handles their relationships for interpolation and differential operators.Source: §3.1.2 · pp. 15–16
3.1.3. Data loading and program execution
- Execution dynamically imports the configured solver, generates the particle object, and repeatedly calls Solver.step() until the simulation ends.The configuration file is supplied to run.py as the program entry point.Source: §3.1.3 · pp. 16–17
3.2. Modules of Major Operators
CraftSPH organizes major SPH operations as reusable modules, letting users combine them to construct solvers with different discretizations and pressure schemes.
3.2.1. Neighbor search
- CraftSPH uses background-grid bucket search and stores source-to-target neighbor relationships and derived quantities for subsequent kernel and differential operations.Separating neighbor generation from later operations makes the processing sequence easier to describe.Source: §3.2.1 · p. 17
3.2.2. Kernel functions
- Users select among seven implemented kernels through configuration, and can add a kernel by defining its value and derivative and registering it.The independent-class interface permits kernel additions without rewriting the solver.Source: §3.2.2 · pp. 17–18
3.2.3. Differential operators
- Differential operators include classical SPH, CSPH, and LSSPH formulations, support matrix forms for implicit schemes, and let new operators use the existing interface.Sparse coefficient matrices reduce memory use as particle counts grow; boundary contributions are stored separately for later Dirichlet or Neumann conditions.Source: §3.2.3 · p. 18
3.2.4. Pressure computation
- Pressure computation supports WCSPH pressure from density variation and implicit linear-system solutions, while users can choose ordinary or symmetric pressure-gradient evaluation.The symmetric pairwise pressure force preserves linear and angular momentum through force antisymmetry.Source: §3.2.4 · p. 19
3.2.5. Stabilization techniques
CraftSPH keeps stabilization and boundary treatments selectable, provides common time-integration updates, and lets users configure output contents.
3.2.5. Stabilization techniques
- Stabilization options include field corrections and particle shifting, called through a common entry point rather than fixed inside the solver.Users can combine techniques according to the problem and selected differential operator, and add new stabilization terms using existing neighbor and kernel information.Source: §3.2.5 · p. 19
3.2.6. Boundary condition treatments
- Boundary handling combines ghost-particle walls, particle-managed inlets and outlets, and free-surface detection, with treatments reusable in explicit and implicit computations.Free-surface classification supports pressure conditions and particle restrictions; boundary data can also enter a coefficient matrix or right-hand-side vector.Source: §3.2.6 · pp. 19–20
3.2.7. Time integration scheme
- The TimeIntegration class provides Euler, second- and fourth-order Runge–Kutta, and Adams–Bashforth state updates, but users implement the full time-advancement sequence.For multistage methods, users must explicitly update particle positions and recompute neighbor searches and kernel functions at each stage.Source: §3.2.7 · pp. 20–21
3.2.8. Data output
- Configuration controls VTK output particle types, variables, names, and intervals, allowing users to select output contents without modifying the solver.Selecting only needed output can reduce data volume, and newly added variables can be included through configuration.Source: §3.2.8 · p. 21
3.3. Solver Construction
A Solver combines common modules and defines the time-evolution procedure, while PyTorch operations support controlled differentiation and SPH–deep-learning integration.
3.3. Solver Construction
- The Solver executes one time step by combining modules in a user-defined procedure that searches neighbors, evaluates forces, integrates particle states, and saves results.Users implement the time-integration procedure, preserving flexibility in solver construction and deep-learning integration.Source: §3.3 · p. 21
- PyTorch automatic differentiation computes loss gradients with respect to selected variables, while discrete operations such as neighbor search may be excluded.Users choose the differentiated range to balance differentiability and computational efficiency.Source: §3.3 · pp. 21–22
- CraftSPH supports deep-learning modules inside SPH solvers and SPH computations inside deep-learning models, with training through losses evaluated through SPH computations.The framework is presented for conventional SPH analysis, differentiable physical simulation, and SPH–deep-learning hybrid models.Source: §3.3 · p. 22
4. EXPERIMENT FOR FORWARD ANALYSIS
Forward tests evaluate CraftSPH’s ISPH-LSSPH solver against analytical solutions, experiments, and benchmark results across internal, multiphase, and free-surface flows. The simulations reproduce key flow behavior across two- and three-dimensional cases, with pressure-history discrepancies and problem-specific accuracy qualifications.
4. EXPERIMENT FOR FORWARD ANALYSIS
- The forward tests used ISPH with LSSPH, demonstrating construction of an advanced solver beyond standard explicit schemes and classical SPH formulations across four fluid-flow benchmarks.The benchmarks cover internal flow, multiphase flow, large-deformation free-surface flow, and three-dimensional flow with complex geometry.Source: §4 · p. 23
4.1. Poiseuille flow
- The predicted Poiseuille velocity profile agreed with the analytical parabolic profile near the walls and remained stable despite continuous inlet addition and outlet removal of particles.Source: §4.1 · p. 24
4.2. Two-dimensional dam break
- Sensor pressure histories reproduced rapid impact rises and subsequent fluctuations, with some instantaneous peak discrepancies; the simulation remained stable for strongly disordered particle distributions.The flow developed local negative-pressure regions, where a similar tendency to suppress tensile instability was observed.Source: §4.2 · pp. 25–26, Figure 4 · p. 26, Figure 5 · p. 27
4.3. Rising bubble
- The Case 1 bubble shape, center-of-mass history, and rise velocity agreed well with the reference, while circularity showed a small discrepancy attributed to particle-based perimeter evaluation.Case 1 was selected because stronger interface deformation in Case 2 would require specialized discretization and stabilization treatments.Source: §4.3 · pp. 27–28, Figure 6 · p. 28
4.4. Three-dimensional dam break
- Pressure histories generally matched experiments after an onset-time discrepancy, while the simulation remained stable and reproduced obstacle impact and complex three-dimensional free-surface motion.The two-dimensional dam-break solver was applied without modifications to the three-dimensional problem.Source: §4.4 · p. 29, Figure 8 · p. 30, Figure 7 · p. 30
5. EXPERIMENT FOR INVERSE ANALYSIS
CraftSPH’s inverse analyses test automatic differentiation in explicit WCSPH and implicit ISPH using high-accuracy LSSPH. Both schemes estimate parameters in two benchmark flows, with accuracy varying by problem and by which solver operations are differentiated.
5. EXPERIMENT FOR INVERSE ANALYSIS
- The inverse-analysis experiments assess parameter estimation from observational data using automatic-differentiation gradients, covering viscosity in Taylor–Green flow and Reynolds number in cavity flow.The study aims to validate gradient computation and applicability to inverse analysis using analytical and established benchmark solutions.Source: §5 · p. 30
- LSSPH is used with both WCSPH and ISPH to test high-accuracy spatial discretization in explicit and implicit formulations, but not every solver operation is differentiated.Neighbor search is excluded from both computational graphs, and the ISPH pressure-Poisson iterative solution is also excluded to avoid potential costs or gradient explosion or vanishing.Source: §5 · pp. 30–31
5.1. Estimation of kinematic viscosity in the Taylor-Green vortex
- In Taylor–Green flow, the viscosity target is estimated by minimizing the mismatch between sensor-interpolated numerical velocities and the analytical solution.Nine sensors are distributed over the domain; each optimization step simulates to 1.0 s, updates viscosity, then restarts from 0 s.Source: §5.1 · p. 31
- WCSPH estimated Taylor–Green kinematic viscosity with 0.19% error, versus 12.13% for ISPH, whose excluded pressure-Poisson solve may reduce gradient accuracy.Both schemes converged to constant estimates and produced final velocity fields generally consistent with theory; total computational times were 14.5 s for WCSPH and 7.7 s for ISPH.Source: §5.1 · p. 32, Table 1 · p. 33
5.2. Estimation of the Reynolds number in lid-driven cavity flow
- For cavity flow, Reynolds number is estimated from benchmark velocity data at Re = 3200 by updating lid velocity during the simulation.The cavity length and viscosity are fixed, parameter updates occur every 0.01 s, and both schemes use the same optimization-step count over the same physical time.Source: §5.2 · pp. 33–34
- Lid-driven cavity Reynolds-number errors were approximately 0.47% for WCSPH and 2.03% for ISPH, with ISPH also sufficiently accurate in this case.The benchmark data came from an independent numerical method; ISPH had lower total computational time, and the authors attribute its sufficient accuracy partly to the small effect of the omitted pressure-Poisson gradient contribution for this optimization target.Source: §5.2 · pp. 34–35, Table 2 · p. 35
6. CONCLUSION
CraftSPH combines modular, differentiable SPH components and supports high-accuracy discretization across varied fluid-flow problems. Its inverse analyses demonstrate automatic differentiation with explicit and implicit schemes, while the paper identifies scope and computational-graph limitations.
6. CONCLUSION
- Inverse analyses demonstrated automatic differentiation with explicit WCSPH and implicit ISPH, with estimation errors depending on the flow problem and scheme.The reported errors were 0.19% and 12.13% for Taylor–Green viscosity estimation, and 0.47% and 2.03% for cavity Reynolds-number estimation, respectively.Source: §6 · p. 36
- CraftSPH is modular: independent numerical components and a unified particle object can be combined to construct and reuse different SPH solvers.A solver developed for a two-dimensional dam-break case was extended to three dimensions without changing its major computational procedures.Source: §6 · p. 36
- Implemented high-accuracy modules supported results across internal flows, strongly deformed free-surface flows, multiphase flows, and three-dimensional problems.The study reports Poiseuille velocity profiles, rising-bubble behavior, and pressure evolution in two- and three-dimensional dam-break problems.Source: §6 · p. 37
- Validation mainly concerns incompressible flows; higher-order differential operators, tensor-valued formulations, and coupling with other numerical methods remain to be developed and validated.Excluding discrete or iterative procedures from automatic differentiation may also affect parameter-estimation accuracy.Source: §6 · p. 37