Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-23 11:02:36 +00:00
Code Issues Packages Projects Releases Wiki Activity

docs(design): Add GPU physics architecture documentation

Browse Source 
- 4 design documents for GPU implementation (2588 lines total)
- gpu_physics_architecture.md: Physics{Elasticity} GPU-first design
- gpu_cpu_backend_architecture.md: Backend selection and dispatch
- gpu_cpu_migration_guide.md: Migration from old API to new
- gpu_elasticity_refactoring.md: Complete refactoring strategy
- Zero CPU-GPU transfer during solve, matrix-free CG
- Elements store geometry, no mesh dependency
- Breaking changes allowed for GPU performance
This commit is contained in:
Jukka Aho
2025-11-12 01:08:22 +02:00
parent 3208904b73
commit 753ddcc9e5
4 changed files with 2588 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/src/book/design/gpu_cpu_backend_architecture.md
+766
View File
@@ -0,0 +1,766 @@
---
title: "GPU-First, CPU-Second Architecture"
date: 2025-11-10
status: "Design"
---
# Backend Architecture: GPU-First, CPU-Second
**Philosophy:** Write once, run on GPU or CPU. Maximum code reuse.
**Strategy:** Julia package extensions for backend-specific code.
## Design Goals
1. **Same API for GPU and CPU** - User code doesn't change
2. **Automatic fallback** - No CUDA? Use CPU automatically
3. **Maximum code reuse** - Shared physics, separate assembly
4. **GPU-first** - Best performance on GPU
5. **CPU-second** - Multithread fallback, not single-thread
## Architecture Overview
```
JuliaFEM.jl/
├── src/
│ ├── physics/
│ │ ├── physics_abstract.jl # Abstract interface (shared)
│ │ ├── physics_elasticity.jl # Physics{Elasticity} (shared)
│ │ └── physics_solvers.jl # CG, Newton (backend-agnostic)
│ └── backends/
│ └── backend_interface.jl # AbstractBackend, CPU(), GPU()
│
├── ext/
│ └── JuliaFEMCUDAExt/
│ ├── cuda_assembly.jl # GPU kernels for assembly
│ ├── cuda_solvers.jl # GPU CG solver
│ └── JuliaFEMCUDAExt.jl # Extension entry point
│
└── Project.toml
[weakdeps]
CUDA = "..."
[extensions]
JuliaFEMCUDAExt = "CUDA"
```
## Package Extensions (Julia 1.9+)
**How it works:**
```toml
# Project.toml
[deps]
LinearAlgebra = "..."
SparseArrays = "..."
[weakdeps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
[extensions]
JuliaFEMCUDAExt = "CUDA"
```
**Behavior:**
- **User has CUDA:** Extension loads automatically → GPU backend available
- **No CUDA:** Extension not loaded → CPU backend only
- **Seamless:** User code same either way
## Abstract Interface (Shared Code)
### Core Types
```julia
# src/physics/physics_abstract.jl
abstract type AbstractBackend end
struct CPU <: AbstractBackend
nthreads::Int
end
CPU() = CPU(Threads.nthreads())
# GPU defined in extension (only if CUDA available)
# struct GPU <: AbstractBackend
# device::CuDevice
# end
abstract type AbstractPhysics{P} end
# Physics{Elasticity} implementation
mutable struct Physics{P} <: AbstractPhysics{P}
name::String
dimension::Int
properties::P
# Geometry and material (backend-agnostic)
body_elements::Vector{Element}
# BCs (backend-agnostic representation)
bc_dirichlet::DirichletBC
bc_neumann::NeumannBC
# Backend-specific data (union over concrete types)
backend_data::Union{Nothing, AbstractBackendData}
end
```
### Boundary Conditions (Shared)
```julia
# src/physics/physics_elasticity.jl
struct DirichletBC
node_ids::Vector{Int}
components::Vector{Vector{Int}}
values::Vector{Vector{Float64}}
end
struct NeumannBC
surface_elements::Vector{Element}
traction::Vector{Vec{3,Float64}}
end
# Add BCs (backend-agnostic)
function add_dirichlet!(
physics::Physics{Elasticity},
node_ids::Vector{Int},
components::Vector{Int},
value::Float64
)
# Same for CPU and GPU
bc = physics.bc_dirichlet
for node in node_ids
push!(bc.node_ids, node)
push!(bc.components, components)
push!(bc.values, fill(value, length(components)))
end
end
```
### Solver Interface (Abstract)
```julia
# src/physics/physics_solvers.jl
"""
Main solver entry point - dispatches to backend
"""
function solve_physics!(
physics::Physics{Elasticity};
backend::AbstractBackend = default_backend(),
time::Float64 = 0.0,
tol = 1e-6,
max_iter = 1000
)
# Initialize backend data
if physics.backend_data === nothing
physics.backend_data = initialize_backend!(physics, backend, time)
end
# Dispatch to backend-specific solver
return solve_backend!(physics, backend, time, tol, max_iter)
end
# Default: Use GPU if available, else CPU
function default_backend()
if @isdefined(CUDA) && CUDA.functional()
return GPU()
else
return CPU()
end
end
```
## CPU Backend (Core - No Extension Needed)
### Assembly with Multithreading
```julia
# src/backends/cpu_assembly.jl
struct CPUBackendData <: AbstractBackendData
# Global arrays (CPU)
nodes::Matrix{Float64} # 3 × n_nodes
elements::Matrix{Int} # 4 × n_elements
E::Vector{Float64}
ν::Vector{Float64}
# BC arrays
is_fixed::Vector{Bool}
prescribed::Vector{Float64}
# Surface loads
surface_nodes::Matrix{Int}
surface_traction::Matrix{Float64}
# Working arrays
f_ext::Vector{Float64}
u::Vector{Float64}
# Threading data
element_stresses::Vector{Vector{SymmetricTensor{2,3,Float64,6}}} # Per thread
end
"""
Compute residual using multithreading
"""
function compute_residual_cpu!(
data::CPUBackendData,
u::Vector{Float64}
)
n_threads = Threads.nthreads()
n_elements = size(data.elements, 2)
# Phase 1: Compute stresses (parallel over elements)
stresses = data.element_stresses
Threads.@threads for elem_idx in 1:n_elements
tid = Threads.threadid()
# Extract element data
conn = data.elements[:, elem_idx]
X = data.nodes[:, conn]
u_elem = [u[3*n-2:3*n] for n in conn]
E = data.E[elem_idx]
ν = data.ν[elem_idx]
# Compute stress (same as GPU kernel, but CPU)
σ = compute_element_stress(X, u_elem, E, ν)
stresses[tid][elem_idx] = σ
end
# Phase 2: Nodal assembly (parallel over nodes)
n_nodes = size(data.nodes, 2)
f_int = zeros(3 * n_nodes)
# Node-to-elements map needed (built during initialization)
Threads.@threads for node_idx in 1:n_nodes
f_node = zero(Vec{3,Float64})
# Loop over elements touching this node
for elem_idx in node_to_elements[node_idx]
# Accumulate force from this element
f_node += compute_nodal_force(node_idx, elem_idx, stresses, data)
end
f_int[3*node_idx-2:3*node_idx] = [f_node[1], f_node[2], f_node[3]]
end
# Residual
r = f_int - data.f_ext
# Apply Dirichlet BC
r[data.is_fixed] .= 0.0
return r
end
```
### CPU CG Solver
```julia
# src/backends/cpu_solvers.jl
function cg_solve_cpu!(
data::CPUBackendData;
tol = 1e-6,
max_iter = 1000
)
n_dofs = length(data.f_ext)
u = data.u
b = data.f_ext
# Initial residual
r = b - compute_residual_cpu!(data, u)
r[data.is_fixed] .= 0.0
p = copy(r)
r_dot_r = dot(r, r)
for iter in 1:max_iter
# Matrix-free: Ap = K*p
Ap = compute_residual_cpu!(data, p)
Ap[data.is_fixed] .= 0.0
alpha = r_dot_r / dot(p, Ap)
u .+= alpha .* p
r .-= alpha .* Ap
r_dot_r_new = dot(r, r)
if sqrt(r_dot_r_new) < tol
@info " CPU CG converged in $iter iterations"
return iter, sqrt(r_dot_r_new)
end
beta = r_dot_r_new / r_dot_r
p .= r .+ beta .* p
r_dot_r = r_dot_r_new
end
@warn "CPU CG did not converge"
return max_iter, sqrt(r_dot_r)
end
```
## GPU Backend (Extension)
### Extension Structure
```julia
# ext/JuliaFEMCUDAExt/JuliaFEMCUDAExt.jl
module JuliaFEMCUDAExt
using JuliaFEM
using CUDA
using Tensors
# Define GPU backend type
struct GPU <: JuliaFEM.AbstractBackend
device::CuDevice
end
GPU() = GPU(CuDevice(0))
# Export to make available when extension loads
export GPU
# GPU-specific data
struct GPUBackendData <: JuliaFEM.AbstractBackendData
nodes::CuArray{Float64,2}
elements::CuArray{Int32,2}
E::CuArray{Float64,1}
ν::CuArray{Float64,1}
is_fixed::CuArray{Bool,1}
prescribed::CuArray{Float64,1}
surface_nodes::CuArray{Int32,2}
surface_traction::CuArray{Float64,2}
f_ext::CuArray{Float64,1}
u::CuArray{Float64,1}
# ... node-to-elements map
end
# Include GPU assembly and solvers
include("cuda_assembly.jl")
include("cuda_solvers.jl")
# Register backend
function JuliaFEM.initialize_backend!(
physics::Physics{Elasticity},
backend::GPU,
time::Float64
)
# Transfer CPU data to GPU
return initialize_gpu_data(physics, time)
end
function JuliaFEM.solve_backend!(
physics::Physics{Elasticity},
backend::GPU,
time, tol, max_iter
)
gpu_data = physics.backend_data::GPUBackendData
# Compute external forces
compute_external_forces_gpu!(gpu_data)
# Solve
iterations, residual = cg_solve_gpu!(gpu_data; tol, max_iter)
# Return CPU array
u_cpu = Array(gpu_data.u)
return (u=u_cpu, iterations=iterations, residual=residual)
end
end # module
```
### GPU Assembly (Same Kernels)
```julia
# ext/JuliaFEMCUDAExt/cuda_assembly.jl
# Same CUDA kernels as before
function compute_element_stresses_kernel!(...)
# ... exact same as gpu_physics_elasticity.jl
end
function nodal_assembly_kernel!(...)
# ... exact same as gpu_physics_elasticity.jl
end
function apply_surface_traction_kernel!(...)
# ... exact same as gpu_physics_elasticity.jl
end
function apply_dirichlet_kernel!(...)
# ... exact same as gpu_physics_elasticity.jl
end
function compute_residual_gpu!(gpu_data::GPUBackendData, u::CuArray)
# ... same implementation
end
```
## User API (Unified)
### Example: Cantilever Beam
```julia
using JuliaFEM
# Optional: Load CUDA support
using CUDA # If installed, GPU backend available
# Create elements (same for CPU and GPU)
body_elements = Element[]
for conn in connectivity
el = Element(Tet4, conn)
update!(el, "geometry", 0.0 => X)
update!(el, "youngs modulus", 0.0 => 210e9)
update!(el, "poissons ratio", 0.0 => 0.3)
push!(body_elements, el)
end
# Create physics (same for CPU and GPU)
physics = Physics(Elasticity, "cantilever", 3)
add_elements!(physics, body_elements)
# Add BCs (same for CPU and GPU)
add_dirichlet!(physics, [1, 2, 3], [1,2,3], 0.0)
surf = Element(Tri3, [4, 5, 6])
update!(surf, "geometry", 0.0 => X_surf)
add_neumann!(physics, surf, Vec{3}((0.0, 0.0, -1e6)))
# Solve - automatic backend selection
result = solve_physics!(physics) # Uses GPU if available, else CPU
# Or explicit backend
result = solve_physics!(physics, backend=GPU()) # Force GPU
result = solve_physics!(physics, backend=CPU(8)) # Force CPU with 8 threads
```
### Backend Selection
```julia
# Automatic (default)
result = solve_physics!(physics)
# Manual
if CUDA.functional()
result = solve_physics!(physics, backend=GPU())
println("Solved on GPU")
else
result = solve_physics!(physics, backend=CPU())
println("Solved on CPU with $(Threads.nthreads()) threads")
end
```
## Code Reuse Analysis
### Shared (80% of code)
**✅ Physics definition:**
- `Physics{Elasticity}` type
- `add_elements!`
- `add_dirichlet!`, `add_neumann!`
- BC data structures
**✅ Solver logic:**
- CG algorithm (same structure)
- Convergence checks
- Newton iteration loop (future)
**✅ Material models:**
- Hooke's law
- Plasticity algorithms
- Stress update
**✅ Integration:**
- Gauss points
- Shape functions
- Jacobian computation
### Backend-Specific (20% of code)
**GPU (in extension):**
- CUDA kernel launches (`@cuda`)
- CuArray operations
- GPU memory management
**CPU (in core):**
- `Threads.@threads` loops
- Regular Array operations
- Thread-local storage
## Performance Characteristics
### GPU Backend
**Strengths:**
- Massive parallelism (1000+ threads)
- Matrix-free (no memory for K)
- Fast for large problems (100K+ DOFs)
**Weaknesses:**
- Small problems: kernel launch overhead
- CPU-GPU transfer (if not careful)
**Best for:** > 10K DOFs
### CPU Backend
**Strengths:**
- No CPU-GPU transfer
- Good single-element performance
- Easier debugging
**Weaknesses:**
- Limited parallelism (typ. 8-16 threads)
- Needs more memory (thread-local arrays)
**Best for:** < 10K DOFs, or no GPU available
### Crossover Point (Projected)
| DOFs | GPU Time | CPU Time (8 threads) | Winner |
|------|----------|----------------------|--------|
| 1K | 50 ms | 20 ms | CPU |
| 10K | 200 ms | 200 ms | TIE |
| 100K | 2 s | 20 s | GPU |
| 1M | 20 s | 200 s | GPU |
## Implementation Plan
### Phase 1: Refactor Current Code (Day 1)
**Goal:** Extract shared interface from `gpu_physics_elasticity.jl`
**Tasks:**
1. Create `src/physics/` directory
2. Move `Physics{Elasticity}` to `physics_elasticity.jl` (no CUDA imports)
3. Move BC structs to shared file
4. Define `AbstractBackend`, `AbstractBackendData`
### Phase 2: CPU Backend (Day 2)
**Goal:** Multithreaded CPU implementation
**Tasks:**
1. Create `src/backends/cpu_assembly.jl`
2. Implement `compute_residual_cpu!` with `Threads.@threads`
3. Implement `cg_solve_cpu!`
4. Test on cantilever beam
### Phase 3: Extension Setup (Day 3)
**Goal:** Move GPU code to extension
**Tasks:**
1. Create `ext/JuliaFEMCUDAExt/` directory
2. Update `Project.toml` with weakdeps and extensions
3. Move GPU kernels to extension
4. Implement `initialize_backend!` and `solve_backend!`
5. Test with and without CUDA
### Phase 4: Unified API (Day 4)
**Goal:** Single entry point for both backends
**Tasks:**
1. Implement `solve_physics!` dispatch
2. Automatic backend detection
3. Update demos to use unified API
4. Document backend selection
### Phase 5: Optimization (Week 2)
**Goal:** Performance tuning
**Tasks:**
1. Benchmark CPU vs GPU
2. Optimize threading (CPU)
3. Optimize kernel parameters (GPU)
4. Add preconditioning (both backends)
## Testing Strategy
### Unit Tests
```julia
@testset "Physics Interface" begin
physics = Physics(Elasticity, "test", 3)
# ... test add_elements!, add_dirichlet!, etc.
end
@testset "CPU Backend" begin
physics = setup_test_problem()
result = solve_physics!(physics, backend=CPU(1))
@test result.residual < 1e-6
end
@testset "GPU Backend" begin
if CUDA.functional()
physics = setup_test_problem()
result = solve_physics!(physics, backend=GPU())
@test result.residual < 1e-6
else
@test_skip "CUDA not available"
end
end
```
### Integration Tests
```julia
@testset "Cantilever CPU vs GPU" begin
physics = setup_cantilever()
result_cpu = solve_physics!(physics, backend=CPU())
if CUDA.functional()
result_gpu = solve_physics!(physics, backend=GPU())
# Solutions should match
@test isapprox(result_cpu.u, result_gpu.u, rtol=1e-6)
end
end
```
### Performance Tests
```julia
using BenchmarkTools
physics = setup_large_problem(n_nodes=100_000)
# CPU
@btime solve_physics!($physics, backend=CPU())
# GPU
if CUDA.functional()
@btime solve_physics!($physics, backend=GPU())
end
```
## Migration Path
### From Current Code
**Old (GPU only):**
```julia
include("src/gpu_physics_elasticity.jl")
using .GPUElasticityPhysics
physics = Physics(Elasticity, "body", 3)
result = solve_elasticity_gpu!(physics)
```
**New (GPU or CPU):**
```julia
using JuliaFEM
physics = Physics(Elasticity, "body", 3)
result = solve_physics!(physics) # Automatic backend
```
**Breaking changes:**
- `solve_elasticity_gpu!` → `solve_physics!`
- Must specify backend explicitly if want GPU-only: `backend=GPU()`
## Benefits
**For Users:**
✅ **Automatic fallback** - No CUDA? No problem, uses CPU
✅ **Same API** - Write once, run anywhere
✅ **Explicit control** - Can force backend if desired
✅ **No dependencies** - CUDA optional, not required
**For Developers:**
✅ **Code reuse** - 80% shared between backends
✅ **Maintainability** - One physics implementation
✅ **Testability** - Can test CPU without GPU
✅ **Extensibility** - Easy to add new backends (Metal, ROCm, etc.)
## Future Backends
**Potential additions:**
- **Metal GPU** (Apple Silicon) - Same extension pattern
- **ROCm** (AMD GPU) - Same extension pattern
- **Distributed CPU** (MPI) - Use same Physics type
- **Intel oneAPI** (Intel GPUs) - Same extension pattern
**All use same `AbstractBackend` interface!**
## File Structure (Final)
```
JuliaFEM.jl/
├── src/
│ ├── JuliaFEM.jl # Main module
│ ├── physics/
│ │ ├── physics_abstract.jl # AbstractPhysics, AbstractBackend
│ │ ├── physics_elasticity.jl # Physics{Elasticity} (shared)
│ │ ├── physics_solvers.jl # solve_physics! dispatch
│ │ └── physics_bcs.jl # BC types (shared)
│ └── backends/
│ ├── backend_interface.jl # AbstractBackendData, CPU type
│ ├── cpu_assembly.jl # CPU multithreaded assembly
│ └── cpu_solvers.jl # CPU CG solver
│
├── ext/
│ └── JuliaFEMCUDAExt/
│ ├── JuliaFEMCUDAExt.jl # Extension entry, GPU type
│ ├── cuda_assembly.jl # GPU kernels
│ └── cuda_solvers.jl # GPU CG solver
│
├── demos/
│ ├── cantilever_unified.jl # Works on CPU or GPU
│ └── backend_comparison.jl # Benchmark both
│
├── test/
│ ├── test_physics_interface.jl # Backend-agnostic tests
│ ├── test_cpu_backend.jl # CPU-specific tests
│ └── test_gpu_backend.jl # GPU-specific tests (skip if no CUDA)
│
└── Project.toml
[deps]
LinearAlgebra = "..."
Tensors = "..."
[weakdeps]
CUDA = "..."
[extensions]
JuliaFEMCUDAExt = "CUDA"
```
## Summary
**Key Insights:**
1. **GPU-first, CPU-second** - Best performance on GPU, reliable fallback on CPU
2. **Package extensions** - CUDA optional, loads automatically if available
3. **Maximum code reuse** - 80% shared, 20% backend-specific
4. **Same user API** - Write once, run on GPU or CPU
5. **Explicit control** - Can force backend if needed
**Next Steps:**
1. Refactor current GPU code to extract shared interface
2. Implement CPU backend with multithreading
3. Set up package extension for GPU
4. Update demos to use unified API
5. Benchmark and optimize
**Result:** Users get best of both worlds - GPU performance when available, CPU fallback always works.
docs/src/book/design/gpu_cpu_migration_guide.md
+774
View File
@@ -0,0 +1,774 @@
---
title: "Migration Guide: GPU-Only to GPU-First/CPU-Second"
date: 2025-11-10
status: "Implementation Plan"
---
# Migration Guide: Adding CPU Backend
**Goal:** Refactor `src/gpu_physics_elasticity.jl` to support both GPU and CPU with maximum code reuse.
**Strategy:** Extract shared interface, implement CPU backend, move GPU to extension.
## Step-by-Step Migration
### Step 1: Extract Shared Types (Day 1, Morning)
**Create:** `src/physics/physics_types.jl`
**Extract from `gpu_physics_elasticity.jl`:**
```julia
# Shared types (no CUDA dependency)
abstract type AbstractBackend end
abstract type AbstractBackendData end
struct CPU <: AbstractBackend
nthreads::Int
end
CPU() = CPU(Threads.nthreads())
# Elasticity properties (no CUDA)
mutable struct Elasticity <: FieldProblem
formulation::Symbol
finite_strain::Bool
geometric_stiffness::Bool
end
# BC types (no CUDA)
struct DirichletBC
node_ids::Vector{Int}
components::Vector{Vector{Int}}
values::Vector{Vector{Float64}}
end
struct NeumannBC
surface_elements::Vector{Element}
traction::Vector{Vec{3,Float64}}
end
# Physics type (backend-agnostic)
mutable struct Physics{P} <: AbstractPhysics{P}
name::String
dimension::Int
properties::P
body_elements::Vector{Element}
bc_dirichlet::DirichletBC
bc_neumann::NeumannBC
backend_data::Union{Nothing, AbstractBackendData}
end
```
**No changes to:**
- `add_elements!`
- `add_dirichlet!`
- `add_neumann!`
These work for both CPU and GPU!
### Step 2: Create Backend Interface (Day 1, Afternoon)
**Create:** `src/physics/physics_interface.jl`
```julia
"""
Initialize backend data from Physics
"""
function initialize_backend! end
"""
Solve physics on specific backend
"""
function solve_backend! end
"""
Main entry point - dispatches to backend
"""
function solve_physics!(
physics::Physics{Elasticity};
backend::AbstractBackend = auto_backend(),
time::Float64 = 0.0,
tol = 1e-6,
max_iter = 1000
)
# Initialize if needed
if physics.backend_data === nothing
physics.backend_data = initialize_backend!(physics, backend, time)
end
# Dispatch to backend
return solve_backend!(physics, backend, time, tol, max_iter)
end
"""
Auto-detect best backend
"""
function auto_backend()
# Check if GPU extension loaded
if isdefined(Main, :CUDA) && Main.CUDA.functional()
return GPU() # Defined in extension
else
return CPU()
end
end
```
### Step 3: Implement CPU Backend (Day 2)
**Create:** `src/backends/cpu_backend.jl`
#### 3a. CPU Data Structure
```julia
struct CPUBackendData <: AbstractBackendData
# Geometry
nodes::Matrix{Float64} # 3 × n_nodes
elements::Matrix{Int32} # 4 × n_elements
n_nodes::Int
n_elements::Int
# Material
E::Vector{Float64}
ν::Vector{Float64}
# BCs
is_fixed::Vector{Bool}
prescribed::Vector{Float64}
# Surface loads
surface_nodes::Matrix{Int32}
surface_traction::Matrix{Float64}
# Node-to-elements (CSR)
node_to_elem_ptr::Vector{Int32}
node_to_elem_data::Vector{Int32}
# Working arrays
f_ext::Vector{Float64}
u::Vector{Float64}
end
```
#### 3b. Initialize CPU Backend
```julia
function initialize_backend!(
physics::Physics{Elasticity},
backend::CPU,
time::Float64
)
@info "Initializing CPU backend ($(backend.nthreads) threads)..."
# Extract data from elements (same as GPU version)
n_elements = length(physics.body_elements)
# Build node map
node_set = Set{Int}()
for el in physics.body_elements
for node in get_connectivity(el)
push!(node_set, node)
end
end
node_list = sort(collect(node_set))
node_map = Dict(node => i for (i, node) in enumerate(node_list))
n_nodes = length(node_list)
# Extract coordinates
nodes = zeros(3, n_nodes)
for el in physics.body_elements
X = el("geometry", time)
conn = get_connectivity(el)
for (local_idx, global_node) in enumerate(conn)
renumbered = node_map[global_node]
nodes[:, renumbered] = X[:, local_idx]
end
end
# Build connectivity
elements = zeros(Int32, 4, n_elements)
for (i, el) in enumerate(physics.body_elements)
conn = get_connectivity(el)
for (j, node) in enumerate(conn)
elements[j, i] = node_map[node]
end
end
# Extract materials
E = [el("youngs modulus", time) for el in physics.body_elements]
ν = [el("poissons ratio", time) for el in physics.body_elements]
# Build BC arrays (same as GPU)
n_dofs = 3 * n_nodes
is_fixed = fill(false, n_dofs)
prescribed = zeros(n_dofs)
bc = physics.bc_dirichlet
for (i, node_id) in enumerate(bc.node_ids)
if haskey(node_map, node_id)
renumbered = node_map[node_id]
for (comp_idx, comp) in enumerate(bc.components[i])
dof = 3 * (renumbered - 1) + comp
is_fixed[dof] = true
prescribed[dof] = bc.values[i][comp_idx]
end
end
end
# Surface loads (same as GPU)
n_surface = length(physics.bc_neumann.surface_elements)
surface_nodes = zeros(Int32, 3, n_surface)
surface_traction = zeros(Float64, 3, n_surface)
for (i, surf_el) in enumerate(physics.bc_neumann.surface_elements)
conn = get_connectivity(surf_el)
for (j, node) in enumerate(conn)
surface_nodes[j, i] = node_map[node]
end
t = physics.bc_neumann.traction[i]
surface_traction[:, i] = [t[1], t[2], t[3]]
end
# Node-to-elements map
node_to_elems = [Int32[] for _ in 1:n_nodes]
for (el_idx, el) in enumerate(physics.body_elements)
for node in get_connectivity(el)
renumbered = node_map[node]
push!(node_to_elems[renumbered], el_idx)
end
end
ptr = zeros(Int32, n_nodes + 1)
ptr[1] = 1
for i in 1:n_nodes
ptr[i+1] = ptr[i] + length(node_to_elems[i])
end
data = vcat(node_to_elems...)
# Create backend data
return CPUBackendData(
nodes, elements, n_nodes, n_elements,
E, ν,
is_fixed, prescribed,
surface_nodes, surface_traction,
ptr, data,
zeros(n_dofs), zeros(n_dofs)
)
end
```
#### 3c. CPU Assembly Functions
```julia
"""
Compute stress at one Gauss point (pure function, no threading)
"""
function compute_stress_at_gp(
X::NTuple{4,Vec{3}},
u::NTuple{4,Vec{3}},
E::Float64,
ν::Float64
)
# Shape derivatives (Tet4)
dN_dxi = (
Vec{3}((-1.0, -1.0, -1.0)),
Vec{3}((1.0, 0.0, 0.0)),
Vec{3}((0.0, 1.0, 0.0)),
Vec{3}((0.0, 0.0, 1.0))
)
# Jacobian
J = dN_dxi[1] ⊗ X[1] + dN_dxi[2] ⊗ X[2] +
dN_dxi[3] ⊗ X[3] + dN_dxi[4] ⊗ X[4]
invJ = inv(J)
# Physical derivatives
dN_dx = (invJ ⋅ dN_dxi[1], invJ ⋅ dN_dxi[2],
invJ ⋅ dN_dxi[3], invJ ⋅ dN_dxi[4])
# Strain
ε = symmetric(dN_dx[1] ⊗ u[1] + dN_dx[2] ⊗ u[2] +
dN_dx[3] ⊗ u[3] + dN_dx[4] ⊗ u[4])
# Stress (Hooke)
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
I = one(ε)
σ = λ * tr(ε) * I + 2μ * ε
return σ, det(J)
end
"""
Compute residual with multithreading
"""
function compute_residual_cpu!(
data::CPUBackendData,
u::Vector{Float64}
)
n_nodes = data.n_nodes
n_elements = data.n_elements
# Phase 1: Internal forces (nodal assembly, parallel over nodes)
f_int = zeros(3 * n_nodes)
Threads.@threads for node_idx in 1:n_nodes
f_node = zero(Vec{3,Float64})
gauss_weight = 1.0 / 24.0 # Tet4
# Loop over elements touching this node
elem_start = data.node_to_elem_ptr[node_idx]
elem_end = data.node_to_elem_ptr[node_idx+1] - 1
for elem_offset in elem_start:elem_end
elem_idx = data.node_to_elem_data[elem_offset]
# Get element data
conn = data.elements[:, elem_idx]
local_node = findfirst(==(node_idx), conn)
# Coordinates
X = tuple([Vec{3}((data.nodes[1, n], data.nodes[2, n], data.nodes[3, n]))
for n in conn]...)
# Displacements
u_elem = tuple([Vec{3}((u[3*n-2], u[3*n-1], u[3*n]))
for n in conn]...)
# Material
E = data.E[elem_idx]
ν = data.ν[elem_idx]
# Compute stress and det(J)
σ, detJ = compute_stress_at_gp(X, u_elem, E, ν)
# Shape derivative for this node
dN_dxi = [Vec{3}((-1.0, -1.0, -1.0)), Vec{3}((1.0, 0.0, 0.0)),
Vec{3}((0.0, 1.0, 0.0)), Vec{3}((0.0, 0.0, 1.0))][local_node]
J = dN_dxi ⊗ X[1] + ... # Recompute (or cache)
invJ = inv(J)
dN_dx = invJ ⋅ dN_dxi
# Accumulate force (4 Gauss points, but Tet4 constant stress)
f_node += (dN_dx ⋅ σ) * (4 * gauss_weight * detJ)
end
f_int[3*node_idx-2] = f_node[1]
f_int[3*node_idx-1] = f_node[2]
f_int[3*node_idx] = f_node[3]
end
# Residual
r = f_int - data.f_ext
# Apply Dirichlet
r[data.is_fixed] .= 0.0
return r
end
```
#### 3d. CPU CG Solver
```julia
function cg_solve_cpu!(
data::CPUBackendData;
tol = 1e-6,
max_iter = 1000
)
u = data.u
b = data.f_ext
# Initial residual
r = b - compute_residual_cpu!(data, u)
r[data.is_fixed] .= 0.0
p = copy(r)
r_dot_r = dot(r, r)
for iter in 1:max_iter
Ap = compute_residual_cpu!(data, p)
Ap[data.is_fixed] .= 0.0
alpha = r_dot_r / dot(p, Ap)
u .+= alpha .* p
r .-= alpha .* Ap
r_dot_r_new = dot(r, r)
if sqrt(r_dot_r_new) < tol
@info " CPU CG converged in $iter iterations"
return iter, sqrt(r_dot_r_new)
end
beta = r_dot_r_new / r_dot_r
p .= r .+ beta .* p
r_dot_r = r_dot_r_new
end
@warn "CPU CG did not converge"
return max_iter, sqrt(r_dot_r)
end
```
#### 3e. CPU Solver Entry Point
```julia
function solve_backend!(
physics::Physics{Elasticity},
backend::CPU,
time::Float64,
tol,
max_iter
)
data = physics.backend_data::CPUBackendData
# Compute external forces (Neumann BC)
@info "Computing external forces (CPU)..."
fill!(data.f_ext, 0.0)
n_surface = size(data.surface_nodes, 2)
for surf_idx in 1:n_surface
# Get nodes
n1, n2, n3 = data.surface_nodes[:, surf_idx]
# Coordinates
X1 = Vec{3}((data.nodes[1, n1], data.nodes[2, n1], data.nodes[3, n1]))
X2 = Vec{3}((data.nodes[1, n2], data.nodes[2, n2], data.nodes[3, n2]))
X3 = Vec{3}((data.nodes[1, n3], data.nodes[2, n3], data.nodes[3, n3]))
# Area
area = 0.5 * norm((X2 - X1) × (X3 - X1))
# Traction
t = Vec{3}((data.surface_traction[1, surf_idx],
data.surface_traction[2, surf_idx],
data.surface_traction[3, surf_idx]))
# Distribute to nodes
force = (area / 3.0) * t
data.f_ext[3*n1-2:3*n1] .+= [force[1], force[2], force[3]]
data.f_ext[3*n2-2:3*n2] .+= [force[1], force[2], force[3]]
data.f_ext[3*n3-2:3*n3] .+= [force[1], force[2], force[3]]
end
# Solve
@info "Solving with CPU CG ($(backend.nthreads) threads)..."
fill!(data.u, 0.0)
iterations, residual = cg_solve_cpu!(data; tol, max_iter)
return (u=copy(data.u), iterations=iterations, residual=residual)
end
```
### Step 4: Move GPU to Extension (Day 3)
**Create:** `ext/JuliaFEMCUDAExt/`
#### 4a. Extension Entry Point
**File:** `ext/JuliaFEMCUDAExt/JuliaFEMCUDAExt.jl`
```julia
module JuliaFEMCUDAExt
using JuliaFEM
using CUDA
using Tensors
# Define GPU backend
struct GPU <: JuliaFEM.AbstractBackend
device::CuDevice
end
GPU() = GPU(CUDA.device())
export GPU
# Include GPU-specific code
include("gpu_backend.jl")
include("gpu_kernels.jl")
end
```
#### 4b. GPU Backend Data
**File:** `ext/JuliaFEMCUDAExt/gpu_backend.jl`
```julia
struct GPUBackendData <: JuliaFEM.AbstractBackendData
# Same fields as CPUBackendData, but CuArrays
nodes::CuArray{Float64,2}
elements::CuArray{Int32,2}
# ... (copy from current gpu_physics_elasticity.jl)
end
# initialize_backend! for GPU
# (copy from current code, rename GPUElasticityData → GPUBackendData)
# solve_backend! for GPU
# (copy from current code)
```
#### 4c. GPU Kernels
**File:** `ext/JuliaFEMCUDAExt/gpu_kernels.jl`
```julia
# Copy all CUDA kernels from current gpu_physics_elasticity.jl:
# - compute_element_stresses_kernel!
# - nodal_assembly_kernel!
# - apply_surface_traction_kernel!
# - apply_dirichlet_kernel!
# - compute_residual_gpu!
# - cg_solve_gpu!
```
### Step 5: Update Project.toml (Day 3)
```toml
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
[weakdeps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
[extensions]
JuliaFEMCUDAExt = "CUDA"
[compat]
julia = "1.9"
CUDA = "5"
```
### Step 6: Update Demos (Day 4)
**Create:** `demos/cantilever_unified.jl`
```julia
using JuliaFEM
# Try to load CUDA
try
using CUDA
println("CUDA available: $(CUDA.functional())")
catch
println("CUDA not available, using CPU")
end
# Create physics (same for CPU or GPU)
physics = Physics(Elasticity, "cantilever", 3)
# ... add elements, BCs
# Solve - automatic backend
result = solve_physics!(physics)
println("Solved with automatic backend")
# Or explicit
if @isdefined(CUDA) && CUDA.functional()
result_gpu = solve_physics!(physics, backend=GPU())
println("GPU result: $(result_gpu.iterations) iterations")
end
result_cpu = solve_physics!(physics, backend=CPU())
println("CPU result: $(result_cpu.iterations) iterations")
# Compare
if @isdefined(result_gpu)
diff = norm(result_gpu.u - result_cpu.u)
println("GPU vs CPU difference: $diff")
end
```
## File Reorganization
### Before (Current)
```
src/
├── JuliaFEM.jl
└── gpu_physics_elasticity.jl (716 lines, all GPU)
```
### After (Target)
```
src/
├── JuliaFEM.jl
├── physics/
│ ├── physics_types.jl # Shared types (100 lines)
│ ├── physics_interface.jl # solve_physics! (50 lines)
│ └── physics_elasticity.jl # add_elements!, add_dirichlet! (100 lines)
└── backends/
└── cpu_backend.jl # CPU implementation (400 lines)
ext/
└── JuliaFEMCUDAExt/
├── JuliaFEMCUDAExt.jl # Extension entry (20 lines)
├── gpu_backend.jl # GPU data structures (100 lines)
└── gpu_kernels.jl # CUDA kernels (400 lines)
```
**Line count:**
- **Shared:** 250 lines (types + interface + elasticity)
- **CPU:** 400 lines (backend + assembly + CG)
- **GPU:** 520 lines (extension + backend + kernels)
- **Total:** 1170 lines (vs 716 before, but now supports CPU!)
## Testing Plan
### Phase 1: Test Extraction (Day 4)
Test that shared types work:
```julia
@testset "Shared Types" begin
physics = Physics(Elasticity, "test", 3)
@test physics.dimension == 3
add_dirichlet!(physics, [1], [1,2,3], 0.0)
@test length(physics.bc_dirichlet.node_ids) == 1
end
```
### Phase 2: Test CPU Backend (Day 4)
```julia
@testset "CPU Backend" begin
physics = setup_simple_problem()
result = solve_physics!(physics, backend=CPU(1))
@test result.residual < 1e-6
end
@testset "CPU Multithreading" begin
physics = setup_simple_problem()
result = solve_physics!(physics, backend=CPU(4))
@test result.residual < 1e-6
@test result.iterations < 500
end
```
### Phase 3: Test GPU Extension (Day 5)
```julia
@testset "GPU Backend" begin
if @isdefined(CUDA) && CUDA.functional()
physics = setup_simple_problem()
result = solve_physics!(physics, backend=GPU())
@test result.residual < 1e-6
else
@test_skip "CUDA not available"
end
end
@testset "CPU vs GPU" begin
if @isdefined(CUDA) && CUDA.functional()
physics = setup_simple_problem()
result_cpu = solve_physics!(physics, backend=CPU())
result_gpu = solve_physics!(physics, backend=GPU())
@test isapprox(result_cpu.u, result_gpu.u, rtol=1e-6)
end
end
```
## Checklist
### Day 1: Extract Shared Code
- [ ] Create `src/physics/physics_types.jl`
- [ ] Move `Elasticity`, `DirichletBC`, `NeumannBC`, `Physics{T}`
- [ ] Create `AbstractBackend`, `CPU` type
- [ ] Test that types work without CUDA
### Day 2: Implement CPU Backend
- [ ] Create `src/backends/cpu_backend.jl`
- [ ] Implement `CPUBackendData`
- [ ] Implement `initialize_backend!(physics, CPU(), time)`
- [ ] Implement `compute_residual_cpu!` with threading
- [ ] Implement `cg_solve_cpu!`
- [ ] Implement `solve_backend!(physics, CPU(), ...)`
- [ ] Test CPU solve on simple problem
### Day 3: Create Extension
- [ ] Create `ext/JuliaFEMCUDAExt/` directory
- [ ] Update `Project.toml` with weakdeps
- [ ] Create `JuliaFEMCUDAExt.jl` entry point
- [ ] Define `GPU` backend type
- [ ] Move GPU kernels to `gpu_kernels.jl`
- [ ] Implement `initialize_backend!(physics, GPU(), time)`
- [ ] Implement `solve_backend!(physics, GPU(), ...)`
- [ ] Test GPU solve with extension
### Day 4: Integration
- [ ] Create `solve_physics!` dispatch
- [ ] Implement `auto_backend()` selection
- [ ] Update `src/JuliaFEM.jl` to include new files
- [ ] Remove old `gpu_physics_elasticity.jl` (or deprecate)
- [ ] Create `demos/cantilever_unified.jl`
- [ ] Test automatic backend selection
### Day 5: Documentation & Testing
- [ ] Update user manual with new API
- [ ] Document backend selection
- [ ] Write CPU backend tests
- [ ] Write GPU extension tests
- [ ] Write comparison tests (CPU vs GPU)
- [ ] Benchmark both backends
## Expected Results
### User Experience
**Without CUDA:**
```julia
julia> using JuliaFEM
julia> physics = Physics(Elasticity, "body", 3)
julia> result = solve_physics!(physics)
[ Info: Initializing CPU backend (8 threads)...
[ Info: CPU CG converged in 430 iterations
```
**With CUDA:**
```julia
julia> using JuliaFEM, CUDA
julia> physics = Physics(Elasticity, "body", 3)
julia> result = solve_physics!(physics)
[ Info: Initializing GPU backend...
[ Info: GPU CG converged in 430 iterations
```
### Performance (Projected)
| Problem Size | CPU (8 threads) | GPU | Speedup |
|--------------|-----------------|-----|---------|
| 1K DOFs | 50 ms | 100 ms | 0.5× |
| 10K DOFs | 500 ms | 300 ms | 1.7× |
| 100K DOFs | 50 s | 3 s | 17× |
| 1M DOFs | 500 s | 20 s | 25× |
## Summary
**Migration strategy:**
1. ✅ **Extract** shared types (Physics, BCs) - no CUDA dependency
2. ✅ **Implement** CPU backend with multithreading
3. ✅ **Move** GPU code to package extension
4. ✅ **Unify** API with `solve_physics!(physics, backend)`
5. ✅ **Test** both backends, compare results
**Result:**
- **Same API** for CPU and GPU
- **Automatic fallback** if no CUDA
- **Maximum code reuse** (250 lines shared)
- **GPU-first** (best performance)
- **CPU-second** (reliable fallback)
**Time estimate:** 5 days full-time work
docs/src/book/design/gpu_elasticity_refactoring.md
+452
View File
@@ -0,0 +1,452 @@
---
title: "GPU Elasticity Refactoring Plan"
date: 2025-11-10
author: "Jukka Aho"
status: "Active"
tags: ["design", "gpu", "elasticity", "refactoring"]
---
## Motivation
Current `src/gpu_elasticity.jl` (477 lines) doesn't follow JuliaFEM architecture:
❌ **Current Issues:**
- Custom `ElasticityPhysics` struct instead of `Problem{Elasticity}`
- Depends on `GmshMesh` (Gmsh-specific)
- BCs defined as vectors, not as `Problem{Dirichlet}` / `Problem{Neumann}`
- Elements don't store coordinates via `update!(element, "geometry", nodes)`
- Solver is monolithic function, not modular assembly
✅ **Target Architecture (JuliaFEM Convention):**
```julia
# Field problem
physics_elasticity = Problem(Elasticity, "body", 3)
add_elements!(physics_elasticity, body_elements)
# Boundary conditions
bc_fixed = Problem(Dirichlet, "fixed end", 3, "displacement")
add_elements!(bc_fixed, boundary_elements)
bc_pressure = Problem(Neumann, "pressure load", 3, "displacement")
add_elements!(bc_pressure, surface_elements)
# Solve
solver = Solver(GPU)
add_problems!(solver, [physics_elasticity, bc_fixed, bc_pressure])
solve!(solver)
```
## Architecture Comparison
### Old (current gpu_elasticity.jl)
```julia
struct ElasticityPhysics
mesh::GmshMesh # ❌ Gmsh-specific
material::ElasticMaterial # ❌ Not Element property
fixed_nodes::Vector{Int} # ❌ Should be Problem{Dirichlet}
pressure_nodes::Vector{Int} # ❌ Should be Problem{Neumann}
pressure_value::Float64
end
# Usage
physics = ElasticityPhysics(mesh, material, fixed_nodes, pressure_nodes, 1e6)
result = solve_elasticity_gpu(physics) # Monolithic
```
### New (JuliaFEM convention)
```julia
# 1. Create elements with geometry (immutable API)
body = Element(Tet4, Lagrange{Tet4, 1}, [1,2,3,4];
fields=(geometry = nodes,
youngs_modulus = 210e9,
poissons_ratio = 0.3))
# 2. Create field problem
physics = Problem(Elasticity, "body", 3)
add_elements!(physics, [body])
# 3. Create boundary conditions (separate problems!)
bc_fixed = Problem(Dirichlet, "fixed", 3, "displacement")
fixed_element = Element(Tri3, Lagrange{Tri3, 1}, [1,2,3];
fields=(geometry = nodes,
displacement_1 = 0.0,
displacement_2 = 0.0,
displacement_3 = 0.0))
add_elements!(bc_fixed, [fixed_element])
bc_pressure = Problem(Neumann, "pressure", 3, "displacement")
pressure_element = Element(Tri3, Lagrange{Tri3, 1}, [4,5,6];
fields=(geometry = nodes,
displacement_traction_force_3 = -1e6))
add_elements!(bc_pressure, [pressure_element])
# 4. Solve with GPU
solver = Solver(GPU, Linear)
add_problems!(solver, [physics, bc_fixed, bc_pressure])
solve!(solver, time)
```
## Key Design Principles
### 1. Elements Store Coordinates
**JuliaFEM Convention:**
```julia
# Element creation with immutable fields
el = Element(Tet4, Lagrange{Tet4, 1}, [1, 2, 3, 4];
fields=(geometry = nodes,
youngs_modulus = 210e9,
poissons_ratio = 0.3))
# Access during assembly
X = el.fields.geometry # Returns 3×4 matrix
E = el.fields.youngs_modulus
```
**Why:** Immutable, type-stable (40-130× faster), GPU-compatible.
### 2. Problem{T} for Physics Types
**JuliaFEM Hierarchy:**
```julia
AbstractProblem
├─ FieldProblem (volume elements)
│ ├─ Elasticity
│ ├─ Heat
│ └─ Truss
└─ BoundaryProblem (surface elements)
├─ Dirichlet
├─ Neumann (implicit, via "traction force" field)
└─ Mortar (contact)
```
**Pattern:**
```julia
mutable struct Elasticity <: FieldProblem
formulation::Symbol # :plane_stress, :plane_strain, :continuum
finite_strain::Bool
geometric_stiffness::Bool
store_fields::Vector{Symbol}
end
# Problem wraps physics type
mutable struct Problem{P<:AbstractProblem}
name::String
dimension::Int # DOFs per node
parent_field_name::String # For BCs: "displacement"
elements::Vector{Element}
dofmap::Dict{Element, Vector{Int}}
assembly::Assembly # K, f, C, g matrices
fields::Dict{String, AbstractField}
properties::P # Elasticity instance
end
```
### 3. Boundary Conditions as Problems
**Dirichlet (fixed displacement):**
```julia
bc = Problem(Dirichlet, "fixed end", 3, "displacement")
# Surface element on boundary (immutable)
fixed_surf = Element(Tri3, Lagrange{Tri3, 1}, [1, 2, 3];
fields=(geometry = nodes,
displacement_1 = 0.0,
displacement_2 = 0.0,
displacement_3 = 0.0))
add_elements!(bc, [fixed_surf])
```
**Neumann (traction/pressure):**
```julia
# Surface element with traction (immutable)
pressure_surf = Element(Tri3, Lagrange{Tri3, 1}, [4, 5, 6];
fields=(geometry = nodes,
displacement_traction_force_3 = -1e6))
# Add to physics problem directly
add_elements!(physics, [pressure_surf])
```
### 4. Material Properties on Elements
```julia
# Material defined per element (immutable - set at construction)
steel_el = Element(Tet4, Lagrange{Tet4, 1}, conn;
fields=(geometry = X,
youngs_modulus = 210e9,
poissons_ratio = 0.3))
aluminum_el = Element(Tet4, Lagrange{Tet4, 1}, conn;
fields=(geometry = X,
youngs_modulus = 69e9,
poissons_ratio = 0.33))
# Heterogeneous materials: create different elements with different properties
```
## Implementation Plan
### Phase 1: Element-Based Storage ✅ (Current Session)
**Goal:** Remove `GmshMesh` dependency, use Element fields with immutable API.
**Changes:**
```julia
# OLD
struct ElasticityPhysics
mesh::GmshMesh # ❌
...
end
# NEW - Immutable elements with fields at construction
body = Element(Tet4, Lagrange{Tet4, 1}, [1,2,3,4];
fields=(geometry = X_nodes,
youngs_modulus = E,
poissons_ratio = ν))
```
**File:** `src/gpu_physics_elasticity.jl` (new immutable-based implementation)
**Tasks:**
1. Remove `GmshMesh`, `ElasticityPhysics`, `ElasticMaterial` structs
2. Accept `Vector{Element}` instead of mesh
3. Extract coordinates from `element.fields.geometry`
4. Extract material from `element.fields.youngs_modulus`
5. Update demo to use immutable element creation
### Phase 2: Problem{Elasticity} Integration
**Goal:** Replace custom solver with `Problem{Elasticity}` pattern.
**Changes:**
```julia
# Create problem
physics = Problem(Elasticity, "continuum", 3)
physics.properties.formulation = :continuum
physics.properties.finite_strain = false
# Add elements
for el in body_elements
add_elements!(physics, el)
end
# Assembly (GPU kernel)
assembly = Assembly()
assemble!(assembly, physics, physics.elements, time)
```
**File:** Extend `src/problems_elasticity.jl` with GPU dispatch
**Tasks:**
1. Add `assemble!(::Assembly, ::Problem{Elasticity}, ::Vector{Element}, time, ::Val{:gpu})`
2. Reuse GPU kernels from Phase 1
3. Keep CPU version untouched (backward compatibility)
4. Support both CPU and GPU via dispatch
### Phase 3: Boundary Conditions as Problems
**Goal:** Replace node vectors with `Problem{Dirichlet}` using immutable elements.
**Changes:**
```julia
# OLD
fixed_nodes = [1, 5, 12, ...]
pressure_nodes = [3, 7, 9, ...]
# NEW - Immutable boundary elements
bc_fixed = Problem(Dirichlet, "fixed end", 3, "displacement")
bc_pressure = Problem(Neumann, "pressure", 3, "displacement")
# Add boundary elements
for node in fixed_node_ids
bc_el = Element(Poi1, Lagrange{Poi1, 1}, [node];
fields=(geometry = nodes[:, node:node],
displacement_1 = 0.0,
displacement_2 = 0.0,
displacement_3 = 0.0))
add_elements!(bc_fixed, bc_el)
end
```
**File:** `src/problems_dirichlet.jl` (GPU dispatch)
**Tasks:**
1. Add GPU assembly for `Problem{Dirichlet}`
2. Integrate with GPU solver
3. Neumann BC: Add `Problem{Neumann}` type (future)
### Phase 4: Solver Integration
**Goal:** Unified GPU/CPU solver interface.
**Changes:**
```julia
# Create solver
solver = Solver(GPU, Linear)
add_problems!(solver, physics, bc_fixed)
# Solve
solver(0.0) # Assemble at time=0
u = solver.assembly.u # Solution
```
**File:** `src/solvers_gpu.jl` (new)
**Tasks:**
1. Create `Solver` struct with GPU field
2. Dispatch `assemble!()` to GPU kernels
3. Matrix-free CG on GPU
4. Export results to `assembly.u`
## Migration Path (Backward Compatibility)
### Keep Old API Working
```julia
# OLD API (deprecated but functional)
physics = ElasticityPhysics(mesh, material, fixed, pressure, 1e6)
result = solve_elasticity_gpu(physics)
# NEW API (preferred)
physics = Problem(Elasticity, "body", 3)
bc = Problem(Dirichlet, "fixed", 3, "displacement")
solver = Solver(GPU, Linear)
solve!(solver)
```
**Strategy:**
1. Keep `gpu_elasticity.jl` as `gpu_elasticity_legacy.jl`
2. Create new `gpu_elasticity_v2.jl` with Problem pattern
3. Add deprecation warnings to old API
4. Update demos to new pattern
5. Remove legacy after 2-3 versions
## Benefits of Refactoring
### Architectural
✅ **Follows JuliaFEM conventions** - Consistent with CPU solver
✅ **Mesh-agnostic** - Works with any mesh format
✅ **Modular** - Field + Boundary problems separate
✅ **Extensible** - Easy to add Neumann, Mortar, etc.
### Practical
✅ **Heterogeneous materials** - Per-element properties
✅ **Multiple BC types** - Dirichlet, Neumann, point loads
✅ **Reusable components** - GPU kernels work for nonlinear too
✅ **Testing** - Use existing test infrastructure
### Performance
✅ **Zero allocation** - Element fields are pre-allocated
✅ **Type-stable** - Material access via Element, not Dict
✅ **GPU-friendly** - Nodal assembly pattern unchanged
## Example: Cantilever Beam (New API)
```julia
using JuliaFEM
# 1. Create mesh (any format)
nodes = [...] # 3×n_nodes matrix
connectivity = [...] # 4×n_elements (Tet4)
# 2. Create body elements (immutable API)
body_elements = Element[]
for i in 1:n_elements
conn = connectivity[:, i]
X = nodes[:, conn]
# Immutable: all fields at construction
el = Element(Tet4, Lagrange{Tet4, 1}, conn;
fields=(geometry = X,
youngs_modulus = 210e9,
poissons_ratio = 0.3))
push!(body_elements, el)
end
# 3. Create field problem
physics = Problem(Elasticity, "cantilever", 3)
physics.properties.formulation = :continuum
add_elements!(physics, body_elements)
# 4. Fixed end (Dirichlet BC - immutable elements)
bc_fixed = Problem(Dirichlet, "fixed end", 3, "displacement")
fixed_node_ids = [1, 2, 5, 9, ...] # X=0 plane
for node_id in fixed_node_ids
X_node = nodes[:, node_id:node_id]
bc_el = Element(Poi1, Lagrange{Poi1, 1}, [node_id];
fields=(geometry = X_node,
displacement_1 = 0.0,
displacement_2 = 0.0,
displacement_3 = 0.0))
add_elements!(bc_fixed, bc_el)
end
# 5. Pressure load (surface traction - immutable)
pressure_element_ids = [...] # Top surface Tri3 elements
for el_id in pressure_element_ids
conn = surface_connectivity[:, el_id]
X = nodes[:, conn]
el = Element(Tri3, Lagrange{Tri3, 1}, conn;
fields=(geometry = X,
displacement_traction_force_3 = -1e6))
add_elements!(physics, el)
end
# 6. Solve with GPU
solver = Solver(GPU, Linear)
add_problems!(solver, physics, bc_fixed)
solver(0.0) # Assemble and solve at time=0
# 7. Results
u = solver.assembly.u
println("Max displacement: ", maximum(abs.(u)))
```
## Status
- **Phase 1:** ✅ DESIGN COMPLETE (this document)
- **Phase 2:** 🔄 IN PROGRESS (implementing element-based storage)
- **Phase 3:** ⏳ PENDING (Problem{Elasticity} integration)
- **Phase 4:** ⏳ PENDING (BC as Problems)
- **Phase 5:** ⏳ PENDING (Solver integration)
## Next Steps
1. **Create `src/gpu_physics_elasticity.jl`** with immutable element API (no GmshMesh)
2. **Update demos** to use immutable element creation with fields at construction
3. **Test** that GPU kernels work with `element.fields.geometry` access
4. **Integrate** with `Problem{Elasticity}` pattern
5. **Add** `Problem{Dirichlet}` GPU assembly
6. **Document** new immutable API in user manual
## References
- `src/problems_elasticity.jl` - CPU elasticity implementation
- `src/problems_dirichlet.jl` - Dirichlet BC pattern
- `src/assembly/problems.jl` - Problem definition
- `test/test_elasticity_*.jl` - Example usage patterns
docs/src/book/design/gpu_physics_architecture.md
+596
View File
@@ -0,0 +1,596 @@
---
title: "GPU Physics Architecture"
date: 2025-11-10
status: "Implemented"
---
# GPU-First Physics{Elasticity} Architecture
**Status:** ✅ Implemented (November 10, 2025)
## Overview
Pure GPU implementation of elasticity solver with:
- **Physics{Elasticity}** (renamed from Problem{Elasticity})
- **All boundary conditions in GPU kernels** (no CPU fallback)
- **Elements store geometry** via `update!(element, "geometry", nodes)`
- **Zero CPU-GPU transfer during solve** (matrix-free CG)
- **Breaking changes allowed** (old API deprecated for GPU performance)
## Design Principles
### 1. Everything on GPU
**Critical Rule:** All assembly, BC application, and solving happens on GPU device.
✅ **Allowed:**
```julia
# GPU kernels for assembly
@cuda threads=256 blocks=n compute_stresses_kernel!(...)
# BC application in device code
@cuda threads=256 blocks=n apply_dirichlet_kernel!(...)
# Matrix-free CG on GPU
K*u ≈ compute_residual_gpu!(u)
```
❌ **Not Allowed:**
```julia
# CPU loops
for element in elements
# assemble on CPU
end
# CPU-GPU transfers in hot path
u_cpu = Array(u_gpu) # SLOW!
process_on_cpu(u_cpu)
u_gpu = CuArray(u_cpu)
```
### 2. Physics{T} (Not Problem{T})
**Naming Convention:** Use "Physics" to emphasize field equations, not workflow.
```julia
# Field physics (volume elements)
physics = Physics(Elasticity, "body", 3)
# Boundary physics (surface elements) - future
bc = Physics(Dirichlet, "fixed", 3, "displacement")
```
**Why:** Clearer separation between:
- **Physics** = PDEs and constitutive laws
- **Solver** = Linear/nonlinear solution strategy
- **Analysis** = Time integration, load stepping
### 3. Elements Store Geometry
**No mesh object dependency.** Each element is self-contained.
```julia
el = Element(Tet4, [1, 2, 3, 4])
# Store coordinates
X = [0.0 1.0 0.0 0.0;
0.0 0.0 1.0 0.0;
0.0 0.0 0.0 1.0]
update!(el, "geometry", 0.0 => X)
# Store material
update!(el, "youngs modulus", 0.0 => 210e9)
update!(el, "poissons ratio", 0.0 => 0.3)
# Access during assembly
X_elem = el("geometry", time)
E = el("youngs modulus", time)
```
**Benefits:**
- Works with any mesh format (Gmsh, Abaqus, Code Aster, hand-coded)
- Heterogeneous materials (per-element properties)
- No global mesh data structure needed
### 4. BCs Integrated in Physics{Elasticity}
**Not separate problems.** BCs are part of the physics definition.
```julia
physics = Physics(Elasticity, "body", 3)
# Add body elements
add_elements!(physics, body_elements)
# Add Dirichlet BC (fixed displacement)
add_dirichlet!(physics, [1, 2, 3], [1,2,3], 0.0) # Fix nodes 1,2,3 (all DOFs)
# Add Neumann BC (surface traction)
surf_el = Element(Tri3, [4, 5, 6])
update!(surf_el, "geometry", 0.0 => X_surf)
add_neumann!(physics, surf_el, Vec{3}((0.0, 0.0, -1e6))) # Pressure
# Solve (BCs applied in GPU kernels)
result = solve_elasticity_gpu!(physics)
```
**GPU Implementation:**
**Dirichlet:**
```julia
# Flag array on GPU
is_fixed::CuArray{Bool,1} # true if DOF is constrained
# Kernel applies BC
function apply_dirichlet_kernel!(r, is_fixed)
dof = threadIdx().x + (blockIdx().x - 1) * blockDim().x
if is_fixed[dof]
r[dof] = 0.0
end
end
```
**Neumann:**
```julia
# Surface element loop in GPU kernel
function apply_surface_traction_kernel!(
f_ext, surface_nodes, surface_traction, nodes
)
surf_idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
# Get surface element geometry
n1, n2, n3 = surface_nodes[:, surf_idx]
X1, X2, X3 = ...
# Compute area and traction
area = 0.5 * norm((X2-X1) × (X3-X1))
force = area * traction / 3 # Lumped to nodes
# Atomic add (allows parallel writes)
CUDA.@atomic f_ext[3*n1-2] += force[1]
...
end
```
## File Structure
```
src/
├── gpu_physics_elasticity.jl # NEW: Pure GPU Physics{Elasticity}
│ ├── Physics{Elasticity} # Main type
│ ├── DirichletBC # BC data structures
│ ├── NeumannBC
│ ├── GPUElasticityData # GPU arrays
│ ├── CUDA kernels # Device code
│ └── solve_elasticity_gpu! # Main solver
│
├── assembly/problems.jl # OLD: Problem{T} (CPU code)
├── problems_elasticity.jl # OLD: CPU elasticity (to be renamed Physics)
└── JuliaFEM.jl # Main module (conditionally includes GPU)
demos/
└── cantilever_physics_gpu.jl # NEW: Example usage
docs/src/book/design/
└── gpu_elasticity_refactoring.md # Design document
```
## API Reference
### Types
```julia
# Physics type
mutable struct Physics{P}
name::String
dimension::Int
properties::P # Elasticity instance
body_elements::Vector{Element}
bc_dirichlet::DirichletBC
bc_neumann::NeumannBC
gpu_data::Union{Nothing, GPUElasticityData}
end
# Elasticity properties
mutable struct Elasticity <: FieldProblem
formulation::Symbol # :continuum
finite_strain::Bool
geometric_stiffness::Bool
end
```
### Constructor
```julia
physics = Physics(Elasticity, name::String, dimension::Int)
```
**Example:**
```julia
physics = Physics(Elasticity, "cantilever beam", 3)
physics.properties.formulation = :continuum
physics.properties.finite_strain = false
```
### Adding Elements
```julia
add_elements!(physics::Physics{Elasticity}, elements::Vector{Element})
add_elements!(physics::Physics{Elasticity}, element::Element)
```
**Example:**
```julia
el = Element(Tet4, [1, 2, 3, 4])
update!(el, "geometry", 0.0 => X_elem)
update!(el, "youngs modulus", 0.0 => 210e9)
update!(el, "poissons ratio", 0.0 => 0.3)
add_elements!(physics, el)
```
### Adding Boundary Conditions
**Dirichlet (Fixed Displacement):**
```julia
add_dirichlet!(
physics::Physics{Elasticity},
node_ids::Vector{Int},
components::Vector{Int}, # [1,2,3] for x,y,z
value::Float64
)
```
**Example:**
```julia
# Fix all DOFs at nodes 1, 2, 3
add_dirichlet!(physics, [1, 2, 3], [1, 2, 3], 0.0)
# Prescribe Z displacement at node 10
add_dirichlet!(physics, [10], [3], -0.01)
```
**Neumann (Surface Traction):**
```julia
add_neumann!(
physics::Physics{Elasticity},
surface_element::Element,
traction::Vec{3,Float64}
)
```
**Example:**
```julia
# Pressure load (1 MPa in -Z direction)
surf_el = Element(Tri3, [4, 5, 6])
update!(surf_el, "geometry", 0.0 => X_surf)
add_neumann!(physics, surf_el, Vec{3}((0.0, 0.0, -1e6)))
```
### Solving
```julia
result = solve_elasticity_gpu!(
physics::Physics{Elasticity};
time::Float64=0.0,
tol=1e-6,
max_iter=1000
)
```
**Returns:**
```julia
(
u = u_cpu, # Displacement vector (CPU array)
iterations = iter, # CG iterations
residual = res # Final residual norm
)
```
**Example:**
```julia
result = solve_elasticity_gpu!(physics; tol=1e-6, max_iter=500)
println("Converged in $(result.iterations) iterations")
println("Max displacement: $(maximum(abs.(result.u)))")
```
## Complete Example
```julia
using JuliaFEM
using CUDA
using Tensors
# 1. Create elements with geometry
nodes = [
0.0 1.0 0.0 1.0
0.0 0.0 1.0 1.0
0.0 0.0 0.0 0.0
]
body_elements = Element[]
for conn in [[1,2,3,4]]
el = Element(Tet4, conn)
X_elem = nodes[:, conn]
update!(el, "geometry", 0.0 => X_elem)
update!(el, "youngs modulus", 0.0 => 210e9)
update!(el, "poissons ratio", 0.0 => 0.3)
push!(body_elements, el)
end
# 2. Create Physics{Elasticity}
physics = Physics(Elasticity, "body", 3)
add_elements!(physics, body_elements)
# 3. Add BCs
add_dirichlet!(physics, [1, 3], [1,2,3], 0.0) # Fix nodes 1,3
surf = Element(Tri3, [2, 4])
update!(surf, "geometry", 0.0 => nodes[:, [2,4]])
add_neumann!(physics, surf, Vec{3}((0.0, 0.0, -1e6)))
# 4. Solve
result = solve_elasticity_gpu!(physics)
println("Max |u|: $(maximum(abs.(result.u)))")
```
## Migration from Old API
### Old (gpu_elasticity.jl - DEPRECATED)
```julia
# OLD API (Gmsh-dependent)
mesh = read_gmsh_mesh("beam.msh")
material = ElasticMaterial(210e9, 0.3)
fixed_nodes = get_surface_nodes(mesh, "FixedEnd")
pressure_nodes = get_surface_nodes(mesh, "PressureSurface")
physics = ElasticityPhysics(mesh, material, fixed_nodes, pressure_nodes, 1e6)
result = solve_elasticity_gpu(physics)
```
### New (gpu_physics_elasticity.jl - CURRENT)
```julia
# NEW API (mesh-agnostic, pure GPU)
# Create elements (from any source: Gmsh, Abaqus, hand-coded)
body_elements = [...] # Elements with geometry
# Physics
physics = Physics(Elasticity, "body", 3)
add_elements!(physics, body_elements)
# BCs (integrated, not separate vectors)
add_dirichlet!(physics, fixed_nodes, [1,2,3], 0.0)
surf_elements = create_surface_elements(pressure_nodes, nodes)
for surf in surf_elements
add_neumann!(physics, surf, Vec{3}((0.0, 0.0, -1e6)))
end
# Solve
result = solve_elasticity_gpu!(physics)
```
**Key Differences:**
1. **No GmshMesh** - Elements store geometry directly
2. **Physics{Elasticity}** instead of ElasticityPhysics struct
3. **BCs integrated** - Not separate node vectors
4. **Type-safe traction** - Vec{3,Float64} not scalar pressure
5. **Modular** - Can add multiple BC types
## GPU Kernel Architecture
### Phase 1: Compute Stresses at Integration Points
```julia
function compute_element_stresses_kernel!(σ_gp, u, nodes, elements, E, ν)
# One thread per Gauss point
# Computes stress from displacement gradient
# Stores σ in global array
end
```
**Parallelism:** `n_elements × 4` threads (Tet4 has 4 Gauss points)
### Phase 2: Nodal Assembly (Internal Forces)
```julia
function nodal_assembly_kernel!(r, σ_gp, nodes, elements, node_to_elems)
# One thread per node
# Loops over touching elements
# Accumulates f_int = ∫ B^T σ dV
end
```
**Parallelism:** `n_nodes` threads
**Key:** No atomics needed (each node owned by one thread)
### Phase 3: Apply Surface Traction (Neumann BC)
```julia
function apply_surface_traction_kernel!(f_ext, surface_nodes, traction, nodes)
# One thread per surface element
# Computes area × traction
# Atomic add to force vector
end
```
**Parallelism:** `n_surface_elements` threads
**Atomics:** Required (multiple surfaces can share nodes)
### Phase 4: Apply Dirichlet BC
```julia
function apply_dirichlet_kernel!(r, is_fixed)
# One thread per DOF
# if is_fixed[dof]: r[dof] = 0.0
end
```
**Parallelism:** `3 × n_nodes` threads
**Fast:** Just array indexing, no computation
### Phase 5: Conjugate Gradient Solver
```julia
function cg_solve_gpu!(gpu_data; tol, max_iter)
while not converged
# Matrix-free: Ap = K*p ≈ residual(p)
Ap = compute_residual_gpu!(gpu_data, p)
# CG update (all on GPU)
alpha = dot(r, r) / dot(p, Ap)
u .+= alpha .* p
r .-= alpha .* Ap
...
end
end
```
**Key:** All vectors stay on GPU (u, r, p, Ap)
## Performance Characteristics
**Strengths:**
✅ **Zero CPU-GPU transfer** during solve
✅ **Matrix-free** (no memory for K matrix)
✅ **Scalable** to millions of DOFs
✅ **Type-stable** (all CuArrays, no Dicts)
**Limitations:**
⏳ **No preconditioning yet** (430 iterations for test case)
⏳ **Tet4 only** (higher-order elements future work)
⏳ **Linear elasticity only** (Newton-Krylov for nonlinear coming)
**Target Performance:**
- **Current:** 430 CG iterations for cantilever (190 nodes, 434 elements)
- **With Jacobi:** ~50-100 iterations
- **With ILU(0):** ~10-20 iterations (GPU ILU challenging)
## Future Work
### Phase 2: Preconditioning (Week 1-2)
**Goal:** Reduce CG iterations 430 → 10-20
**Approach:** Chebyshev-Jacobi preconditioner (GPU-friendly)
```julia
# Extract diagonal (matrix-free)
function extract_diagonal_kernel!(diag, nodes, elements, E, ν)
# Compute K_ii by finite difference
# Or: Assemble diagonal of B^T D B
end
# Apply preconditioner
M_inv = Diagonal(1 ./ diag)
CG(M_inv * K, M_inv * f)
```
### Phase 3: Nonlinear Elasticity (Week 3-6)
**Goal:** Newton-Krylov framework for plasticity
**Changes:**
1. **Material state at integration points**
```julia
state_old::CuArray{PlasticState,1} # Per Gauss point
state_new::CuArray{PlasticState,1}
```
2. **Newton loop**
```julia
while norm(R) > tol
R = compute_residual_gpu!(u, state_old)
K_tangent = approximate_jacobian_gpu(u, state_old)
Δu = cg_solve_gpu!(K_tangent, -R)
u += Δu
update_state!(state_new, u) # Trial state
end
commit_state!(state_old, state_new) # Accept converged state
```
3. **Line search**
```julia
α = backtracking_line_search_gpu(u, Δu, R)
u += α * Δu
```
### Phase 4: Higher-Order Elements (Month 2-3)
**Goal:** Tri6, Tet10, Quad8, Hex20 support
**Approach:**
- Extend topology module (already has Tet10 definition)
- Add integration point data (more Gauss points)
- Generalize kernels (variable nodes per element)
### Phase 5: Contact Mechanics (Month 4-6)
**Goal:** Mortar contact on GPU
**Challenges:**
- Pairing algorithm (spatial search on GPU)
- Gap function evaluation (surface projections)
- Contact constraints (augmented Lagrangian)
## Testing
**Unit Tests:**
```julia
# Test BC application
@testset "Dirichlet BC" begin
physics = Physics(Elasticity, "test", 3)
add_dirichlet!(physics, [1], [1,2,3], 0.0)
initialize_gpu!(physics)
@test all(physics.gpu_data.is_fixed[[1,2,3]] .== true)
end
```
**Integration Tests:**
```julia
# Test cantilever beam
@testset "Cantilever GPU" begin
physics = setup_cantilever()
result = solve_elasticity_gpu!(physics)
@test result.iterations < 500
@test result.residual < 1e-6
end
```
**Validation Tests:**
```julia
# Compare to analytical solution
@testset "Beam bending" begin
u_fem = solve_cantilever_gpu()
u_analytical = beam_theory(E, I, L, P)
@test isapprox(u_fem[end], u_analytical, rtol=0.1)
end
```
## References
- **Implementation:** `src/gpu_physics_elasticity.jl`
- **Demo:** `demos/cantilever_physics_gpu.jl`
- **Design:** `docs/src/book/design/gpu_elasticity_refactoring.md`
- **Old API:** `src/gpu_elasticity.jl` (DEPRECATED)
## Authors
- Jukka Aho (original author, maintainer)
- Refactored: November 10, 2025
## License
MIT (see LICENSE.md)