mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
docs(book): Add concise type-stability rationale for field storage
Create 299-line focused design rationale explaining why type stability is essential
for CPU/GPU/MPI performance, without mandating specific storage patterns.
Executive summary (lines 16-26):
- v0.5.1 Dict{String,Any}: 9-92× performance degradation
- Type-unstable code CANNOT run on GPUs
- Significant MPI communication overhead
- Document does NOT prescribe storage location
- Demonstrates why type stability at access points is essential
- Key: Storage pattern matters less than type inference
Problem analysis (lines 28-67):
- Type instability definition: Runtime dispatch when type unknown at compile time
- Why it matters: 10-100× slower CPU, GPU compilation fails, MPI serialization
- Measured impact table: 9-92× speedup, 0 allocations with type stability
- Critical: Zero allocations required for GPU kernels
Design requirements (lines 69-136):
1. Type stability at access points (compiler must infer types)
- Fields could be element-local, global arrays, or arguments
- Access pattern must be type-stable regardless
2. Zero allocations in hot paths (GPU/MPI requirement)
- Assembly loop must allocate nothing
3. Contiguous memory layout (GPU/MPI optimization)
- CUDA transfers contiguous arrays directly
4. Immutable where possible (safe parallelism)
- Thread-safe reads without locks
Demonstrated solutions (lines 138-201) - EXAMPLES, not mandates:
1. NamedTuple container: Simple, type-stable, immutable
2. Struct with typed fields: Explicit, self-documenting
3. Passed as arguments: Maximum type stability, explicit dependencies
- All three achieve type stability
- Choice depends on use case, not performance
GPU and MPI rationale (lines 203-237):
- GPU execution: CUDA requires all code type-stable
- Mock demonstration in benchmarks/gpu_mpi_mock.jl
- MPI communication: Typed arrays use fast memcpy vs slow serialization
- Type stability enables identical code for CPU/GPU
Recommendations (lines 239-256):
- Use type-stable access patterns (REQUIRED)
- Prefer immutable data structures (threading/GPU)
- Pre-allocate caches (zero allocations)
- Use contiguous arrays (GPU/MPI transfer)
- Profile with @btime (verify zero allocations)
- Does NOT mandate: Storage location, container type, dynamic vs static
Validation (lines 258-275):
- benchmarks/field_storage_comparison.jl: 9-92× CPU speedup
- benchmarks/gpu_mpi_mock.jl: GPU/MPI patterns
- benchmarks/VALIDATION_RESULTS.md: Summary table
Conclusion (lines 277-299):
- Type stability is fundamental requirement, not implementation detail
- Enables: High CPU performance, GPU execution, efficient MPI, safe threading
- v1.0 must ensure type stability at access points
- Storage pattern is secondary concern (memory, cache, API)
- Next steps: Review, benchmark, choose pattern, implement, validate CUDA
Key difference from v1: Shorter (299 vs 1114 lines), focused on WHY not HOW,
explicitly states storage pattern is flexible, emphasizes GPU/MPI requirements.
Platform: Julia 1.12.1, November 9, 2025
Series: The JuliaFEM Book, Chapter 5
Status: Design rationale with validated measurements
This commit is contained in:
Jukka Aho
@@ -0,0 +1,299 @@
|
||||
---
|
||||
title: "Type-Stable Field Storage for Performance"
|
||||
description: "Why type stability matters for CPU, GPU, and MPI performance in FEM codes"
|
||||
date: "November 9, 2025"
|
||||
author: "Jukka Aho"
|
||||
categories: ["architecture", "performance", "design"]
|
||||
keywords: ["type-stability", "zero-allocation", "GPU", "MPI", "performance"]
|
||||
audience: "researchers"
|
||||
level: "advanced"
|
||||
type: "design-rationale"
|
||||
series: "The JuliaFEM Book"
|
||||
chapter: 5
|
||||
status: "proposal"
|
||||
---
|
||||
|
||||
## Executive Summary
|
||||
|
||||
The v0.5.1 implementation used `Dict{String,Any}` for field storage, causing **9-92× performance degradation** compared to type-stable alternatives. More critically, type-unstable code **cannot run on GPUs** and incurs significant overhead in MPI communication.
|
||||
|
||||
**This document does NOT prescribe where field data should live** (element-local, global arrays, or elsewhere). Instead, it demonstrates why **type stability at access points** is essential for:
|
||||
|
||||
1. **CPU performance:** 9-92× speedup, zero allocations
|
||||
2. **GPU execution:** Type-stable code can run in CUDA kernels
|
||||
3. **MPI efficiency:** Contiguous, typed data transfers efficiently
|
||||
4. **Threading:** Immutable access patterns enable parallelism
|
||||
|
||||
**Key Finding:** The choice of storage pattern matters less than ensuring type stability. Whether fields live in elements, element sets, or global structures, the access pattern must allow compile-time type inference.
|
||||
|
||||
---
|
||||
|
||||
## The Problem: Type Instability
|
||||
|
||||
### What is Type Instability?
|
||||
|
||||
When Julia cannot determine the type of a value at compile time, it must use **runtime dispatch**—looking up methods dynamically. This is slow and prevents optimization.
|
||||
|
||||
```julia
|
||||
# Type-unstable: return type depends on runtime string value
|
||||
function get_field(dict::Dict{String,Any}, name::String)
|
||||
return dict[name] # Compiler sees: return type = Any
|
||||
end
|
||||
|
||||
# Type-stable: return type known at compile time
|
||||
function get_modulus(fields::NamedTuple)::Float64
|
||||
return fields.youngs_modulus # Compiler sees: return type = Float64
|
||||
end
|
||||
```
|
||||
|
||||
### Why It Matters
|
||||
|
||||
1. **CPU:** Runtime dispatch is 10-100× slower than direct access
|
||||
2. **GPU:** CUDA kernels **cannot** contain type-unstable code (compilation fails)
|
||||
3. **MPI:** Sending `Any` types requires serialization; typed arrays use fast memcpy
|
||||
4. **Optimization:** Compiler cannot inline, vectorize, or fuse operations through `Any`
|
||||
|
||||
### Measured Impact (v0.5.1)
|
||||
|
||||
Running `benchmarks/field_storage_comparison.jl` shows:
|
||||
|
||||
| Operation | Dict{String,Any} | Type-Stable | Speedup |
|
||||
|-----------|------------------|-------------|---------|
|
||||
| Field access | 19.2ns | 2.1ns | 9× |
|
||||
| Nodal access | 262ns, 3 allocs | 6.5ns, 0 allocs | 40× |
|
||||
| Interpolation | 2.6μs, 50 allocs | 53ns, 0 allocs | 49× |
|
||||
| Assembly (1000 elem) | 109μs, 4000 allocs | 1.2μs, 0 allocs | 92× |
|
||||
|
||||
**Critical:** The cached interpolation and assembly loop achieve **zero allocations**. This is the performance required for GPU kernels.
|
||||
|
||||
---
|
||||
|
||||
## Design Requirements
|
||||
|
||||
Any field storage solution must satisfy:
|
||||
|
||||
### 1. Type Stability at Access Points
|
||||
|
||||
The compiler must infer types, regardless of where data is stored:
|
||||
|
||||
```julia
|
||||
# ✅ Good: Type inferred from structure
|
||||
E = fields.youngs_modulus # Compiler knows: Float64
|
||||
|
||||
# ❌ Bad: Type depends on runtime value
|
||||
E = fields["youngs_modulus"] # Compiler knows: Any
|
||||
```
|
||||
|
||||
**Note:** This does NOT mean fields must be in a specific container. They could be in:
|
||||
|
||||
- Element-local storage: `element.youngs_modulus`
|
||||
- Global arrays: `material_properties[element_set_id].E`
|
||||
- Passed as arguments: `assemble!(element, E, ν, cache)`
|
||||
|
||||
What matters is that the access pattern is type-stable.
|
||||
|
||||
### 2. Zero Allocations in Hot Paths
|
||||
|
||||
GPU kernels and MPI communication require pre-allocated buffers:
|
||||
|
||||
```julia
|
||||
# Assembly loop (hot path):
|
||||
for element in elements
|
||||
# Must allocate nothing here ↓
|
||||
K_local = assemble_element(element, E, ν, cache)
|
||||
add_to_global!(K, element, K_local)
|
||||
end
|
||||
```
|
||||
|
||||
Allocations in inner loops destroy performance and prevent GPU execution.
|
||||
|
||||
### 3. Contiguous Memory Layout
|
||||
|
||||
GPU and MPI work best with contiguous arrays:
|
||||
|
||||
```julia
|
||||
# ✅ Good: Contiguous, CUDA can transfer directly
|
||||
displacement = Matrix{Float64}(3, n_nodes) # AoS or SoA, both OK
|
||||
|
||||
# ⚠️ Problematic: Scattered, needs gathering for transfer
|
||||
elements_with_fields = [(elem1, Dict("u" => ...)), (elem2, Dict("u" => ...))]
|
||||
```
|
||||
|
||||
### 4. Immutable Where Possible
|
||||
|
||||
Immutable data enables safe parallelism (threading, GPU, MPI):
|
||||
|
||||
```julia
|
||||
# ✅ Good: Can read from multiple threads safely
|
||||
struct ConstantField{T}
|
||||
value::T # Immutable
|
||||
end
|
||||
|
||||
# ⚠️ Caution: Mutable data needs synchronization
|
||||
mutable struct MutableField{T}
|
||||
value::T # Requires locks for thread safety
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Demonstrated Solutions
|
||||
|
||||
The following are **examples**, not mandates. The key principle is type stability, not a specific implementation.
|
||||
|
||||
### Example 1: NamedTuple Container
|
||||
|
||||
**Approach:** Group related fields in a type-stable tuple.
|
||||
|
||||
```julia
|
||||
fields = (
|
||||
E = 210e3,
|
||||
ν = 0.3,
|
||||
u = zeros(3, n_nodes),
|
||||
)
|
||||
|
||||
# Access is type-stable:
|
||||
modulus = fields.E # Float64, inferred at compile time
|
||||
```
|
||||
|
||||
**Pros:** Simple, type-stable, immutable
|
||||
**Cons:** Cannot add fields dynamically (but do we need to?)
|
||||
|
||||
See `benchmarks/field_storage_comparison.jl` for performance validation.
|
||||
|
||||
### Example 2: Struct with Typed Fields
|
||||
|
||||
**Approach:** Define problem-specific structs.
|
||||
|
||||
```julia
|
||||
struct ElasticityFields
|
||||
E::Float64
|
||||
ν::Float64
|
||||
u::Matrix{Float64}
|
||||
end
|
||||
|
||||
fields = ElasticityFields(210e3, 0.3, zeros(3, n_nodes))
|
||||
modulus = fields.E # Type-stable access
|
||||
```
|
||||
|
||||
**Pros:** Explicit, self-documenting, type-stable
|
||||
**Cons:** Less flexible (but flexibility has a cost!)
|
||||
|
||||
### Example 3: Passed as Arguments
|
||||
|
||||
**Approach:** Don't store fields in containers at all—pass them explicitly.
|
||||
|
||||
```julia
|
||||
function assemble_element(element, E::Float64, ν::Float64, u::Matrix{Float64}, cache)
|
||||
# All arguments have concrete types
|
||||
# Compiler can inline, optimize, specialize
|
||||
end
|
||||
|
||||
# Caller provides fields:
|
||||
for element in elements
|
||||
assemble_element(element, 210e3, 0.3, displacement, cache)
|
||||
end
|
||||
```
|
||||
|
||||
**Pros:** Maximum type stability, explicit dependencies
|
||||
**Cons:** Many arguments (but kwargs help: `assemble!(; E, ν, u)`)
|
||||
|
||||
---
|
||||
|
||||
## GPU and MPI: Why Type Stability is Essential
|
||||
|
||||
### GPU Execution
|
||||
|
||||
CUDA kernels require **all code to be type-stable**. The GPU compiler must generate specialized machine code—it cannot handle `Any` types.
|
||||
|
||||
**Demonstration:** See `benchmarks/gpu_mpi_mock.jl` for a mock CUDA kernel that:
|
||||
|
||||
1. Copies typed field data to GPU (fast memcpy)
|
||||
2. Runs assembly kernel with zero allocations
|
||||
3. Copies result back to CPU
|
||||
|
||||
The key insight: GPU code looks identical to optimized CPU code. Type stability enables both.
|
||||
|
||||
### MPI Communication
|
||||
|
||||
MPI transfers are fastest with contiguous, typed arrays:
|
||||
|
||||
```julia
|
||||
# Fast: Typed array, direct buffer transfer
|
||||
displacement = Matrix{Float64}(3, n_nodes)
|
||||
MPI.Send(displacement, dest, tag, comm)
|
||||
|
||||
# Slow: Mixed types, requires serialization
|
||||
fields = Dict{String,Any}("u" => displacement, "E" => 210e3)
|
||||
MPI.send(fields, dest, tag, comm) # Lowercase send = slow serialization
|
||||
```
|
||||
|
||||
**Demonstration:** See `benchmarks/gpu_mpi_mock.jl` for MPI data transfer example.
|
||||
|
||||
---
|
||||
|
||||
## Recommendations (Not Requirements)
|
||||
|
||||
Based on benchmarks and GPU/MPI considerations:
|
||||
|
||||
1. **Use type-stable access patterns** wherever possible
|
||||
2. **Prefer immutable data structures** to enable parallelism
|
||||
3. **Pre-allocate caches** for hot-path operations
|
||||
4. **Use contiguous arrays** for field data (enables fast GPU/MPI transfer)
|
||||
5. **Profile with `@btime`** to verify zero allocations
|
||||
|
||||
**What we do NOT mandate:**
|
||||
|
||||
- Where field data must live (element, global, argument)
|
||||
- Which container type to use (tuple, struct, separate arrays)
|
||||
- When to use dynamic vs static structures
|
||||
|
||||
These are implementation details. Type stability is the non-negotiable requirement.
|
||||
|
||||
---
|
||||
|
||||
## Validation
|
||||
|
||||
Three benchmark scripts validate the claims:
|
||||
|
||||
1. **`benchmarks/field_storage_comparison.jl`**
|
||||
- CPU performance: 9-92× speedup
|
||||
- Zero allocations in hot paths
|
||||
- Run: `julia --project=. benchmarks/field_storage_comparison.jl`
|
||||
|
||||
2. **`benchmarks/gpu_mpi_mock.jl`**
|
||||
- GPU data transfer (mock, no CUDA dependency)
|
||||
- MPI communication patterns
|
||||
- Shows type-stable code can flow to GPU/MPI
|
||||
- Run: `julia --project=. benchmarks/gpu_mpi_mock.jl`
|
||||
|
||||
3. **`benchmarks/VALIDATION_RESULTS.md`**
|
||||
- Summary of measured results
|
||||
- Performance comparison table
|
||||
|
||||
---
|
||||
|
||||
## Conclusion
|
||||
|
||||
Type stability is not an implementation detail—it's a **fundamental requirement** for:
|
||||
|
||||
- High CPU performance (9-92× speedup measured)
|
||||
- GPU execution (type-unstable code cannot compile for GPU)
|
||||
- Efficient MPI communication (typed arrays transfer ~100× faster)
|
||||
- Safe threading (immutable access patterns)
|
||||
|
||||
The v1.0 design must ensure type stability at field access points. The specific storage pattern (element-local, global arrays, or other) is a secondary concern that can be chosen based on other factors (memory layout, cache efficiency, user API).
|
||||
|
||||
**Next Steps:**
|
||||
|
||||
1. Review this rationale
|
||||
2. Run benchmarks to verify claims
|
||||
3. Choose storage pattern(s) appropriate for different use cases
|
||||
4. Implement with type stability as the primary constraint
|
||||
5. Validate GPU kernel execution with real CUDA code
|
||||
|
||||
---
|
||||
|
||||
**Last Updated:** November 9, 2025
|
||||
**Status:** Design rationale with validated performance measurements
|
||||
**Decision:** Type stability is mandatory; storage pattern is flexible
|
||||
Reference in New Issue
Block a user