v1.19.0

Warp v1.19.0

Warp v1.19 is a regular feature release with updates across the library. Highlights include reusable sparse assembly, new sparse preconditioners, Surface Nets meshing, CPU performance improvements, and experimental texture-sampling gradients.

New features

Reuse sparse structure across updates

In time stepping or a nonlinear FEM solve, a matrix's coefficients can change repeatedly while the mesh connectivity and degrees of freedom stay fixed. Reusing the nonzero pattern avoids sorting and merging the same row and column indices at every assembly. warp.sparse.bsr_set_from_triplets(..., build_mapping=True) returns a BsrTripletMap for subsequent value updates with bsr_set_from_mapped_values() (#1923).

import warp as wp
import warp.sparse as sparse

rows = wp.array([0, 0, 1], dtype=wp.int32)
cols = wp.array([0, 0, 1], dtype=wp.int32)
values = wp.array([1.0, 2.0, 3.0], dtype=wp.float32)
A = sparse.bsr_zeros(2, 2, block_type=wp.float32)
# New in 1.19: build a reusable map for later value updates.
mapping = sparse.bsr_set_from_triplets(
    A, rows, cols, values, prune_numerical_zeros=False, build_mapping=True
)
values.assign([4.0, 5.0, 6.0])
sparse.bsr_set_from_mapped_values(A, mapping, values)
print(A.values.numpy()[:A.nnz_sync()])  # [9. 6.]

Build mappings with compact topology and prune_numerical_zeros=False, and keep the destination topology unchanged. The default gather reduction avoids atomics. reduction="atomic" permits an accumulation order that can vary between runs.

For repeated FEM bilinear assembly, select the corresponding construction policy, then reuse the same output matrix:

import warp as wp
import warp.fem as fem

@fem.integrand
def mass(s: fem.Sample, u: fem.Field, v: fem.Field):
    return u(s) * v(s)

mesh = fem.Grid2D(res=wp.vec2i(2, 2))
domain = fem.Cells(mesh)
space = fem.make_polynomial_space(mesh, degree=1)
fields = {"u": fem.make_trial(space, domain=domain), "v": fem.make_test(space, domain=domain)}
# New in 1.19: cache the assembly structure for repeated integrations.
A = fem.integrate(mass, fields=fields, bsr_options={"construction": "mapped"})
fem.integrate(mass, fields=fields, output=A,
              bsr_options={"construction": "mapped", "capacity": "reuse"})
print(A.nrow, A.ncol, A.nnz_sync())  # 9 9 49

construction="mapped" uses an atomic-free reduction. construction="mapped_atomic" avoids temporary triplet-value storage, with a floating-point accumulation order that can vary between runs. Rebuild the cached map with capacity="auto" after changing the assembly setup or modifying the matrix elsewhere. See the reuse requirements for the full constraints. These policies support CUDA graph replay after setup.

Surface Nets meshing

Surface Nets offers a different vertex layout from Marching Cubes that can produce more evenly spaced vertices and fewer skinny triangles. This makes it worth trying for visualization and downstream mesh processing, including pipelines that use only triangles. warp.geometry.IsoSurfaceNets places vertices inside cells by averaging the points where the surface crosses their edges. It also produces quads directly, and shares the IsoSurfaceBase interface with Marching Cubes so you can compare extractors using the same scalar field (#1614).

import numpy as np
import warp as wp
import warp.geometry

x, y, z = np.meshgrid(*([np.linspace(-1.0, 1.0, 5)] * 3), indexing="ij")
field = wp.array(np.sqrt(x*x + y*y + z*z) - 0.6, dtype=wp.float32)
# New in 1.19: extract a quad mesh with Surface Nets.
vertices, indices = warp.geometry.IsoSurfaceNets.extract(
    field, lower=(-1.0, -1.0, -1.0), upper=(1.0, 1.0, 1.0), topology="quad"
)
print(vertices.shape[0], indices.shape[0] // 4)  # 32 30

Inputs must be three-dimensional wp.float32 arrays with at least two nodes on every axis. Gradients propagate from vertex positions to field values, with cell classification and topology held fixed. Surfaces contained inside the grid are closed. Surfaces extending beyond it produce open meshes and may stop short of the grid boundary.

The isosurface example now accepts --method and --topology to select the extractor and face type.

Sparse solver preconditioners

FSAI preconditioning

Symmetric positive-definite systems that are difficult for Jacobi preconditioning now have another option. warp.optim.linear.FSAI builds a factorized sparse approximate inverse on the matrix's device and works as the M argument to conjugate-gradient or conjugate-residual solves. For fixed compact topology, reuse_pattern=True lets update() refit the factors after values change (#2039).

import warp as wp
import warp.sparse as sparse
from warp.optim.linear import FSAI, cg

rows = wp.array([0, 0, 1, 1], dtype=wp.int32)
cols = wp.array([0, 1, 0, 1], dtype=wp.int32)
A = sparse.bsr_from_triplets(2, 2, rows, cols, wp.array([4.0, 1.0, 1.0, 3.0], dtype=wp.float32))
# New in 1.19: build an FSAI preconditioner with reusable factor sparsity.
M = FSAI(A, reuse_pattern=True)
sparse.bsr_scale(A, 2.0)
M.update()
b = wp.array([1.0, 2.0], dtype=wp.float32)
x = wp.zeros_like(b)
cg(A, b, x, M=M, tol=1.0e-6, use_cuda_graph=False)
print(x.numpy().round(6))  # [0.045455 0.318182]

Store both halves of the symmetric matrix. FSAI does not infer a missing triangle or globally check symmetry and positive definiteness. Inputs may use wp.float16, wp.float32, or wp.float64. Construction synchronizes and must run outside graph capture. Updates keep the selected factor pattern, so rebuild if convergence deteriorates. After an eager warm-up update, updates can be captured in CUDA graphs. Inspect update_status after replay for numerical failures. FSAI does not support automatic differentiation.

Block-Jacobi preconditioning

Try block-Jacobi when scalar Jacobi converges slowly and your matrix groups strongly coupled unknowns into blocks. Treating each diagonal block as a small coupled system lets the preconditioner account for relationships that scalar Jacobi ignores. Choose "block_jacobi_direct", "block_jacobi_sequential", or "block_jacobi_tile" through warp.optim.linear.preconditioner() (#1891).

import warp as wp
import warp.sparse as sparse
from warp.optim.linear import cg, preconditioner

blocks = wp.array([wp.mat22(4.0, 1.0, 1.0, 3.0)], dtype=wp.mat22)
A = sparse.bsr_diag(blocks)
# New in 1.19: precondition with full diagonal blocks.
M = preconditioner(A, "block_jacobi_direct")
b = wp.array([1.0, 2.0], dtype=wp.float32)
x = wp.zeros_like(b)
cg(A, b, x, M=M, tol=1.0e-6, use_cuda_graph=False)
print(x.numpy().round(6))  # [0.090909 0.636364]

All three require a square BSR matrix with square blocks. The sequential and tile variants require symmetric positive-definite diagonal blocks. The tile variant supports wp.float32 and wp.float64. Scalar 1×1 blocks fall back to ordinary diagonal preconditioning.

Textures

Texture sampling and gradients

Important

This is an experimental feature. The API may change without a formal deprecation cycle.

Optimization can now differentiate a sampled texture value with respect to its sampling coordinates and mipmap level of detail. Coordinate gradients require linear spatial filtering, and level-of-detail gradients require linear mip filtering. Closest filtering produces zero coordinate gradients, and LOD gradients are zero where the level of detail is clamped. Texture data remains non-differentiable (#1301).

import numpy as np
import warp as wp

@wp.kernel
def sample(tex: wp.Texture1D, u: wp.array[float], out: wp.array[float]):
    out[0] = wp.texture_sample(tex, u[0], dtype=float)

tex = wp.Texture1D(np.array([0.0, 2.0, 4.0, 6.0], dtype=np.float32),
                   normalized_coords=False, filter_mode=wp.TextureFilterMode.LINEAR)
u = wp.array([1.0], dtype=wp.float32, requires_grad=True)
out = wp.zeros(1, requires_grad=True)
with wp.Tape() as tape:
    wp.launch(sample, dim=1, inputs=[tex, u], outputs=[out])
# New in 1.19: differentiate the texture sampling coordinate.
tape.backward(grads={out: wp.ones_like(out)})
print(u.grad.numpy())  # [2.]

Closest sampling is also available again for wp.int32 and wp.uint32 textures, after its rejection in Warp 1.18. These samples convert numerically to wp.float32 without normalization. Integer textures default to closest filtering and reject linear filtering (#2065):

import numpy as np
import warp as wp

@wp.kernel
def sample(tex: wp.Texture1D, out: wp.array[float]):
    out[0] = wp.texture_sample(tex, 0.75, dtype=float)

# Restored in 1.19: closest sampling of 32-bit integer textures.
tex = wp.Texture1D(np.array([10, 123], dtype=np.int32))
out = wp.zeros(1)
wp.launch(sample, dim=1, inputs=[tex], outputs=[out])
print(out.numpy())  # [123.]

Tile programming

Optional TF32 matrix multiplication

Tile matrix multiplication now offers TensorFloat-32 (TF32), which uses lower input precision than wp.float32 while retaining wp.float32 accumulation and storage. Pass precision="tf32" to wp.tile_matmul() (#937). Selecting precision="tf32" may not improve performance with the current implementation.

import warp as wp

@wp.kernel
def multiply(A: wp.array2d[float], B: wp.array2d[float], C: wp.array2d[float]):
    a = wp.tile_load(A, shape=(16, 16))
    b = wp.tile_load(B, shape=(16, 16))
    # New in 1.19: request TF32 input precision.
    wp.tile_store(C, wp.tile_matmul(a, b, precision="tf32"))

A = wp.ones((16, 16), dtype=wp.float32, device="cuda:0")
C = wp.zeros_like(A)
wp.launch_tiled(multiply, dim=1, inputs=[A, A], outputs=[C], block_dim=128, device="cuda:0")
print(C.numpy()[0, :4])  # [16. 16. 16. 16.]

TF32 requires the cuBLASDx path and a CUDA compile target of compute capability 8.0 or higher. CPU execution and unsupported CUDA paths use wp.float32 arithmetic and emit a compilation warning. Adjoint matrix products use the same precision choice. Tensor Core use depends on the workload and tile configuration.

We aim to improve TF32 performance in future releases by optimizing Warp's internal data layouts and how intermediate results are accumulated (#1965).

Python and compiler workflows

Array subclass support

Initialized subclasses of wp.array, wp.indexedarray, wp.fabricarray, and wp.indexedfabricarray can now be used wherever their parent array family is accepted. Your application can define its own array classes and Python helper methods while passing instances directly to kernels, function overloads, wp.map(), and struct fields (#1723).

import warp as wp

class TaggedArray(wp.array):
    def __init__(self, data, *, tag, **kwargs):
        super().__init__(data, **kwargs)
        self.tag = tag

@wp.kernel
def square(x: wp.array[float]):
    i = wp.tid()
    x[i] = x[i] * x[i]

x = TaggedArray([2.0, 3.0], tag="sensor", dtype=wp.float32)
# New in 1.19: pass the initialized subclass as an array argument.
wp.launch(square, dim=2, inputs=[x])
print(x.tag, x.numpy().tolist())  # sensor [4.0, 9.0]

Python metadata and helper methods are unavailable inside kernels. New allocations and views are not guaranteed to preserve the subclass or its metadata. Kernel specialization uses the parent family, dtype, and dimensionality. See the array subclass guide for ownership and descriptor requirements.

Avoid object-array overhead in NumPy

Passing a Warp scalar type to NumPy constructors such as np.zeros() or np.full() previously created an object array, adding allocation and conversion overhead during setup. With NumPy 2.4 or newer, these calls now allocate native numeric arrays, avoiding that slowdown when preparing simulation data (#2038):

import numpy as np
import warp as wp

# New in 1.19: use a Warp scalar dtype directly in NumPy 2.4+.
print(np.zeros(3, dtype=wp.float32).dtype)  # float32

wp.bfloat16 is excluded. With older NumPy versions, use wp.dtype_to_numpy().

Count set bits

Integer bit masks can now be counted with wp.bit_count(), at Python scope or inside kernels. It accepts signed and unsigned scalar integers from 8 to 64 bits and returns wp.int32. Signed inputs count the fixed-width two's-complement representation, including the sign bit (#2071).

import warp as wp

# New in 1.19: count bits in each integer's fixed-width representation.
print(wp.bit_count(wp.uint32(45)), wp.bit_count(wp.int8(-1)))  # 4 8

Performance improvements

CPU execution and compilation

CPU compiler improvements can reduce simulation time without changes to kernels. Warp 1.19 enables loop and component vectorization and uses more aggressive helper inlining (#1954, #1994). Gains depend on the workload and CPU, and cold compilation can take longer.

Independent local checks compared the Warp 1.18.0 release wheel with a source build of Warp 1.19.0 on Linux, using one logical CPU of an AMD EPYC 9B45. These are complete Newton simulation-step timings after warmup, excluding compilation, model setup, and rendering:

Newton CPU example	Warp 1.18 (ms/frame)	Warp 1.19 (ms/frame)	Execution speedup
Hanging cloth (32 × 16, VBD)	189	123	1.53×
100 quadruped worlds (XPBD)	15	10.2	1.47×
Cable Y-junction (VBD)	9.85	4.57	2.16×
Hanging softbody (VBD)	278	114	2.45×

Both versions ran identical Newton examples. Results vary by workload and CPU.

Compilation tradeoff: More aggressive CPU optimization can increase compilation time. Reusing compiled kernels amortizes that cost, while short jobs or workloads that frequently compile new kernels may take longer overall.

Other performance improvements

• Sparse CUDA operations: Matrix construction uses more compact sort keys, and matrix transposition avoids sorting work (#1986, #1971).
• Tile kernels: Kernels that allocate tiles inside loops or @wp.func calls, including some backward passes, spend less time accessing shared tile storage (#2085).
• Python overhead: Repeated declarations using wp.static() reuse parsed source, and Python-scope builtins have lower first-call latency for uncached argument signatures (#1486, #1990).

API Capture updates

Source fallback and loading from memory

Important

This is an experimental feature. The API may change without a formal deprecation cycle.

Saved CUDA graphs can now optionally retain generated source for recompilation when a packaged binary cannot load on the replay GPU. Set source_fallback=True in wp.capture_save() to include that source alongside the compiled kernels. The file-based loader tries compatible binaries before recompiling the retained source (#1837).

This option adds CUDA source JIT compilation on the deployment target when fallback is needed. It does not export a multi-architecture CUDA fatbin.

Source fallback requires the exact Warp version used for export. See the API Capture guide for compatibility limits.

When saving an existing API Capture graph:

# New in 1.19: include source when saving an existing API Capture graph.
wp.capture_save(graph, "demo", source_fallback=True)

For a graph saved as "demo", deploy both the .wrp file and its companion module directory:

demo.wrp            # Recorded operations, snapshots, and module metadata
demo_modules/
    *.cubin / *.ptx  # Compiled CUDA kernels
    *.cu            # Generated source retained with source_fallback=True

For standalone C++ replay, see the loader requirements and example.

Related updates:

• Load graphs embedded in other artifacts. If a capture ships inside a model or application archive, wp.capture_load_from_bytes() loads its .wrp and module buffers without first writing them out as separate files (#1724).
• Initialize mesh-handle arrays directly. wp.array(data, dtype=wp.handle) now accepts Python sequences and NumPy arrays. Use this type for mesh IDs that must be remapped when a saved graph is loaded (#2114).

Given those saved artifacts, load and replay from buffers:

from pathlib import Path
import warp as wp

modules = {p.name: p.read_bytes() for p in Path("demo_modules").iterdir() if p.is_file()}
# New in 1.19: load the saved graph and its modules from memory.
graph = wp.capture_load_from_bytes(Path("demo.wrp").read_bytes(), modules, device="cuda:0")
wp.capture_launch(graph)

Loading from bytes requires compatible packaged binaries. Source recompilation and persistent replay-device caches are available through the file-based loader.

New and updated examples

• The distributed 3-D FDTD example simulates a Luneburg lens across multiple GPUs, using NCCL through nccl4py for halo exchange. It requires mpi4py, an MPI implementation, and the nccl4py extra matching the CUDA major version (#2086).
• The isosurface example compares Marching Cubes and Surface Nets, including quad output (#1614).

Breaking changes

Removals in this release

• The sparse masked argument is removed. Use topology="masked" in bsr_set_from_triplets(), bsr_assign(), bsr_set_transpose(), bsr_axpy(), and bsr_mm(). Remove explicit masked=False to retain compact behavior (#1537).

- sparse.bsr_set_from_triplets(A, rows, cols, values, masked=True)
+ sparse.bsr_set_from_triplets(A, rows, cols, values, topology="masked")
- sparse.bsr_mm(A, B, masked=False)
+ sparse.bsr_mm(A, B)

• Unused Marching Cubes arguments and compatibility attributes are removed. These were retained for compatibility and already had no effect on extraction. Output arrays were sized dynamically, and extraction already used the input field's device. Omit max_verts, max_tris, and device when constructing IsoSurfaceMarchingCubes or its deprecated wp.MarchingCubes alias, and omit max_verts and max_tris from resize() calls. Remove accesses to the corresponding compatibility attributes, as well as id and runtime (#1594).

- extractor = warp.geometry.IsoSurfaceMarchingCubes(nx, ny, nz, max_verts=V, max_tris=T, device=device)
+ extractor = warp.geometry.IsoSurfaceMarchingCubes(nx, ny, nz)

.wrp file format bumped to version 18

API Capture files now use .wrp format version 18 to fix mesh-ID relocation from handle arrays and struct arrays. Recapture graphs saved in earlier formats. For mesh IDs that must remain valid after loading in another process, use wp.handle (#2114):

- ids = wp.array([mesh.id], dtype=wp.uint64)
+ ids = wp.array([mesh.id], dtype=wp.handle)

When writing mesh IDs during capture, use a kernel that receives the mesh as a wp.handle argument. IDs written with array.fill_() during capture are not remapped on load, even in a wp.handle array.

Deterministic kernels reject unsafe access patterns

Deterministic mode now rejects kernels that mix deterministic atomic updates with direct reads, writes, or ordinary atomics on the same array (#1817). For example, this kernel raises wp.WarpDeterminismError when compiled:

import warp as wp

@wp.kernel(module="unique", module_options={"deterministic": wp.DeterministicMode.RUN_TO_RUN})
def add_then_read(values: wp.array[float], observed: wp.array[float]):
    i = wp.tid()
    wp.atomic_add(values, i, 1.0)
    # New in 1.19: reading the deterministic atomic target is rejected.
    observed[i] = values[i]

wp.load_module(module=add_then_read.module, device="cpu")

To make floating-point accumulation repeatable, Warp records the atomic updates and applies them in a fixed order after the kernel finishes. The read of values[i] above would therefore happen before the increment is applied, despite appearing after it in the source. Rejecting the kernel prevents silently changing its result. Use separate input and output arrays, or move the reads into a later kernel launch.

Parallel CUDA backward passes that both modify an array and accumulate into its gradient must also be split. Deterministic compilation failures now report a cause and suggested fix through wp.WarpDeterminismError, a RuntimeError subclass. Update exception handlers accordingly:

- except wp.WarpCodegenError:
+ except wp.WarpDeterminismError:

Array subclass launch and gradient changes

Kernel launches now use the data, dtype, and shape stored by the Warp array itself, including for subclasses. An overridden NumPy or CUDA array interface can no longer substitute different data or a different layout. If your subclass relied on such reinterpretation, pass an explicit view or converted array instead. For example, if a stores float32 data shaped (4, 3) and kernel expects wp.array[wp.vec3] (#1723):

- wp.launch(kernel, dim=4, inputs=[a])
+ wp.launch(kernel, dim=4, inputs=[a.view(wp.vec3)])

Use a value conversion when you intend to change numerical dtype, rather than reinterpret its bytes.

Automatic differentiation now also uses the subclass's .grad array and retain_grad setting, as for a regular wp.array. This can change gradients when kernels overwrite outputs or update arrays in place. Review those cases against the migration guide.

Corrected loop-unrolling limit

An off-by-one fix means stepped range() loops that previously slipped past max_unroll may now compile as dynamic loops. This can affect gradients because dynamic-loop differentiation has limitations.

Migration: Review stepped range() loops in kernels that rely on unrolling. Where needed, raise max_unroll to at least the loop's actual iteration count (#1980).

For example, range(0, 33, 2) has 17 iterations, so preserving unrolling requires a limit of at least 17:

- wp.set_module_options({"max_unroll": 16})
+ wp.set_module_options({"max_unroll": 17})

Custom FEM interpolation API changes

FEM interpolation converts requested physical field values into coefficients stored by a function space. This conversion previously failed to reproduce even constant vector fields in some Nédélec and Raviart-Thomas spaces. The fix changes FunctionSpace.dof_value() to receive the weights and coordinate transform needed to compute those coefficients (#2062):

- dof = space.dof_value(value, node_weight, local_value_map)
+ dof = space.dof_value(value, node_dual_weight, local_inverse_value_map)

Migration for custom implementations: If you implement your own FunctionSpace, provide local_inverse_value_map_inner() and local_inverse_value_map_outer(). Custom vector-valued BasisSpace or ShapeFunction implementations must also provide make_node_dual_weight(). Warp's built-in spaces already provide these methods. The same update also fixes interpolation in quadratic and cubic B-spline spaces.

Announcements

Upcoming removals

• Python 3.10 support is deprecated and will be removed in Warp 1.20. Python 3.10 reached end of life on October 1, 2026. Runtime and source-build warnings now identify this change. Upgrade to Python 3.11 or newer to continue receiving Warp updates (#2128).

• The per-environment sequence form of warp.fem.Nanogrid.from_environment_voxels() and warp.fem.AdaptiveNanogrid.from_environment_voxels() will be removed in Warp 1.20. This form was deprecated in Warp 1.16. Pass flat points, point_envs, and env_count instead, with cell_levels for adaptive grids.

• warp.sparse.BsrMatrix.copy_nnz_async() will be removed in Warp 1.20. Use nnz_sync() to read the exact host-side block count, and notify_nnz_changed() after directly modifying sparse storage metadata.

Windows on ARM64 wheels

Important

Windows on ARM64 wheel support is experimental.

• Prebuilt Windows ARM64 wheels are now published on PyPI. They support CPU execution and CUDA acceleration, and are built without MathDx. Nightly builds do not include these wheels. Use a release wheel to try this platform.

Acknowledgments

We also thank the following contributors from outside the core Warp development team:

• @and2049 for fixing Clang CUDA compilation with kernel launch bounds (#1968).
• @bodapatisaikrishna for fixing BVH construction and cell lookups for two-dimensional deformed FEM geometries (#1993).
• @eigenvivek for adding texture-sampling gradients (#1301).
• @EylonKrause for fixing out-of-bounds histogram merging in large tile sorts (#1748).
• @jaredmejia for fixing nested-loop adjoints and empty strided-loop iteration (#2048, #2046).
• @songyinggoh for fixing tile-map adjoints and static builtin hashing (#2031, #2014).
• @v0ropaev for validating host vector components and fixing multidimensional array casts (#2076, #2123).
• @zhihuidu-amd for fixing CUDA tile reductions with partially filled thread blocks and improving native compilation portability (#1996).
• @AntoineRichard for fixing biased Poisson sampling on CUDA devices and ARM CPUs (#2083).
• @dilenvin for adding block-Jacobi preconditioning for sparse linear solvers (#1891).
• @nvjonwong for adding reusable FEM assembly with fixed topology (#1923).
• @jnskkmhr for fixing multi-environment nanogrid cell lookups (#2017).
• @maxkra15 for improving sparse matrix transposition and first-sample FEM interpolation (#1970, #1971).

---

For a complete list of changes, see the full changelog.