From 753ddcc9e5a39eea5140a8d75f9c31ab341dab1c Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 12 Nov 2025 01:08:22 +0200 Subject: [PATCH] docs(design): Add GPU physics architecture documentation - 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 --- .../design/gpu_cpu_backend_architecture.md | 766 +++++++++++++++++ .../book/design/gpu_cpu_migration_guide.md | 774 ++++++++++++++++++ .../book/design/gpu_elasticity_refactoring.md | 452 ++++++++++ .../book/design/gpu_physics_architecture.md | 596 ++++++++++++++ 4 files changed, 2588 insertions(+) create mode 100644 docs/src/book/design/gpu_cpu_backend_architecture.md create mode 100644 docs/src/book/design/gpu_cpu_migration_guide.md create mode 100644 docs/src/book/design/gpu_elasticity_refactoring.md create mode 100644 docs/src/book/design/gpu_physics_architecture.md diff --git a/docs/src/book/design/gpu_cpu_backend_architecture.md b/docs/src/book/design/gpu_cpu_backend_architecture.md new file mode 100644 index 0000000..849ab2b --- /dev/null +++ b/docs/src/book/design/gpu_cpu_backend_architecture.md @@ -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. diff --git a/docs/src/book/design/gpu_cpu_migration_guide.md b/docs/src/book/design/gpu_cpu_migration_guide.md new file mode 100644 index 0000000..d7dc2c9 --- /dev/null +++ b/docs/src/book/design/gpu_cpu_migration_guide.md @@ -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 diff --git a/docs/src/book/design/gpu_elasticity_refactoring.md b/docs/src/book/design/gpu_elasticity_refactoring.md new file mode 100644 index 0000000..4f1368c --- /dev/null +++ b/docs/src/book/design/gpu_elasticity_refactoring.md @@ -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 diff --git a/docs/src/book/design/gpu_physics_architecture.md b/docs/src/book/design/gpu_physics_architecture.md new file mode 100644 index 0000000..ffaed35 --- /dev/null +++ b/docs/src/book/design/gpu_physics_architecture.md @@ -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)