chore: removed files from demos

Removed files from demos directory:

- demos/assembly_comparison_simple.jl
- demos/cantilever_beam_demo.jl
- demos/cantilever_beam_simple.jl
- demos/cantilever_cpu_comparison.jl
- demos/cantilever_gmsh_gpu.jl
- demos/cantilever_physics_gpu.jl
- demos/gpu_assembly_poc_tensors.jl
- demos/gpu_assembly_poc.jl
- demos/gpu_assembly_tet10.jl
- demos/gpu_elementset_demo.jl
- demos/GPU_KERNEL_PLAN.md
- demos/gpu_mpi_demo.jl
- demos/gpu_mpi_mock.jl
- demos/gpu_nodal_assembly_demo.jl
- demos/gpu_only_demo.jl
- demos/krylov_mpi_gpu_demo.jl
- demos/newton_krylov_anderson_cpu.jl
- demos/nodal_assembly_cpu.jl
- demos/nodal_assembly_gpu.jl
- demos/README_GPU_MPI.md
- demos/README_GPU_POC.md
- demos/README_KRYLOV_DEMO.md
- demos/README_TENSORS_CORRECTION.md
- demos/test_tet10_cpu.jl
This commit is contained in:
Jukka Aho
2025-12-12 22:00:41 +02:00
parent cee336a476
commit 9c5229bf4a
24 changed files with 0 additions and 8116 deletions
-296
View File
@@ -1,296 +0,0 @@
# GPU Kernel Plan: Newton-Krylov-Anderson
## What We Learned from CPU Reference
The complete solver pipeline (see `newton_krylov_anderson_cpu.jl`) has these operations:
### Outer Loop: Newton Iterations
```julia
for iter in 1:max_newton_iter
r = assemble_residual(u) # ← GPU KERNEL 1
du = gmres_solve(J, -r) # ← GMRES loop (see below)
u_new = u + α * du # ← GPU vector op
u = anderson_step(u_new, r) # ← CPU (small vectors)
end
```
### Middle Loop: GMRES Iterations
```julia
for j in 1:max_gmres_iter
w = J * v_j # ← GPU KERNEL 1 (matvec via FD)
# Arnoldi orthogonalization # ← cuBLAS (dot products)
# Givens rotations # ← CPU (small matrices)
end
```
### Inner Operation: Matrix-Free Matvec
```julia
function apply_jacobian(v)
r_pert = assemble_residual(u + ε*v) # ← GPU KERNEL 1
return (r_pert - r) / ε
end
```
## GPU Kernel Requirements
### **KERNEL 1: Element Residual Assembly** (The Workhorse)
**What it does:**
```julia
for elem in elements # ← parallelized over GPU threads
for gp in gauss_points
# 1. Compute Jacobian
J = Σ dN_i X_i
# 2. Physical derivatives
dN_dx = J^-1 · dN_dξ
# 3. Strain
ε = sym(Σ dN_i u_i)
# 4. Stress with plasticity
(σ, state_new) = return_mapping(ε, state_old)
# 5. Nodal forces
f_i = dN_i · σ
# 6. Scatter to global (atomics!)
atomic_add!(r_global[dofs], f_i)
end
end
```
**Inputs:** (all on GPU)
- `nodes::CuArray{Float64, 2}` - Shape (3, n_nodes)
- `elements::CuArray{Int32, 2}` - Shape (4, n_elements) for Tet4
- `u::CuArray{Float64}` - Solution vector
- `states::CuArray{PlasticState}` - Plastic states per GP
**Outputs:** (all on GPU)
- `r::CuArray{Float64}` - Residual vector
- `states::CuArray{PlasticState}` - Updated states (in-place)
**Key challenge:** Manual indexing inside kernel (no Tensors.jl)
### **KERNEL 2: Vector Operations** (Trivial)
These are standard cuBLAS operations:
- `y = α*x + β*y` (axpy)
- `dot(x, y)`
- `norm(x)`
Already provided by CUDA.jl
## The GPU Implementation Strategy
### Phase 1: Single Kernel Test
```julia
# Test just the residual assembly kernel
r_cpu = assemble_residual_cpu(u)
r_gpu = assemble_residual_gpu(u) # ← Implement this!
@test r_cpu r_gpu
```
### Phase 2: Matrix-Free Matvec Test
```julia
# Test Jacobian-vector product
Jv_cpu = apply_jacobian_cpu(u, v)
Jv_gpu = apply_jacobian_gpu(u, v) # ← Just calls kernel twice
@test Jv_cpu Jv_gpu
```
### Phase 3: GMRES on GPU
```julia
# Krylov.jl supports CuArrays!
using Krylov, CUDA
du_gpu = CuVector(zeros(n))
gmres!(du_gpu, J_gpu, -r_gpu) # ← Should work with CuArrays
```
### Phase 4: Complete Newton-Krylov-Anderson on GPU
```julia
u_gpu = CuVector(u)
states_gpu = CuArray(states)
for iter in 1:max_iter
r_gpu = assemble_residual_gpu!(u_gpu, states_gpu) # ← GPU
gmres!(du_gpu, J_gpu, -r_gpu) # ← GPU
u_gpu .+= du_gpu # ← GPU
u_cpu = Vector(u_gpu) # ← Copy to CPU
u_cpu = anderson_step(u_cpu, Vector(r_gpu)) # ← CPU
u_gpu = CuVector(u_cpu) # ← Copy to GPU
end
```
**Note:** Anderson runs on CPU (small vectors, needs least-squares). This is fine!
## Key Implementation Details
### 1. Plastic State on GPU
```julia
# CPU version (Tensors.jl)
struct PlasticState
ε_p::SymmetricTensor{2,3,Float64,6}
α::Float64
end
# GPU version (manual indexing)
struct PlasticStateGPU
ε_p::NTuple{6,Float64} # Voigt notation: (ε11, ε22, ε33, ε12, ε13, ε23)
α::Float64
end
```
### 2. Return Mapping on GPU
Must be **completely manual** - no Tensors.jl:
```julia
function return_mapping_gpu(ε_total, state_old, E, ν, σ_y)
# All operations in Voigt notation
# ε = (ε11, ε22, ε33, ε12, ε13, ε23)
# Trial stress (manual Hooke's law)
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
ε_e = ε_total .- state_old.ε_p
tr_ε = ε_e[1] + ε_e[2] + ε_e[3]
σ_trial = (
λ * tr_ε + 2μ * ε_e[1], # σ11
λ * tr_ε + 2μ * ε_e[2], # σ22
λ * tr_ε + 2μ * ε_e[3], # σ33
2μ * ε_e[4], # σ12
2μ * ε_e[5], # σ13
2μ * ε_e[6] # σ23
)
# Deviatoric stress (manual)
p = (σ_trial[1] + σ_trial[2] + σ_trial[3]) / 3
s = (σ_trial[1] - p, σ_trial[2] - p, σ_trial[3] - p,
σ_trial[4], σ_trial[5], σ_trial[6])
# von Mises stress
σ_eq = sqrt(3/2 * (s[1]^2 + s[2]^2 + s[3]^2 + 2*(s[4]^2 + s[5]^2 + s[6]^2)))
# Yield check
f = σ_eq - σ_y
if f <= 0
return (σ_trial, state_old, false)
else
# Radial return
Δγ = f / (3μ)
β = 1 - 2μ * Δγ / σ_eq
σ_new = (
β * s[1] + p,
β * s[2] + p,
β * s[3] + p,
β * s[4],
β * s[5],
β * s[6]
)
# Update plastic strain
n = s ./ σ_eq
Δε_p = Δγ .* n
ε_p_new = state_old.ε_p .+ Δε_p
α_new = state_old.α + Δγ
state_new = PlasticStateGPU(ε_p_new, α_new)
return (σ_new, state_new, true)
end
end
```
### 3. Element Kernel Structure
```julia
function element_residual_kernel!(r_global, nodes, elements, u, states, E, ν, σ_y)
elem_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if elem_idx <= size(elements, 2)
# Extract element nodes
n1, n2, n3, n4 = elements[:, elem_idx]
X1 = (nodes[1, n1], nodes[2, n1], nodes[3, n1])
X2 = (nodes[1, n2], nodes[2, n2], nodes[3, n2])
X3 = (nodes[1, n3], nodes[2, n3], nodes[3, n3])
X4 = (nodes[1, n4], nodes[2, n4], nodes[3, n4])
# Extract displacements
u1 = (u[3*n1-2], u[3*n1-1], u[3*n1])
u2 = (u[3*n2-2], u[3*n2-1], u[3*n2])
u3 = (u[3*n3-2], u[3*n3-1], u[3*n3])
u4 = (u[3*n4-2], u[3*n4-1], u[3*n4])
# Shape derivatives (constant for Tet4)
dN1 = (-1.0, -1.0, -1.0)
dN2 = (1.0, 0.0, 0.0)
dN3 = (0.0, 1.0, 0.0)
dN4 = (0.0, 0.0, 1.0)
# Jacobian (manual tensor product and sum)
J11 = dN1[1]*X1[1] + dN2[1]*X2[1] + dN3[1]*X3[1] + dN4[1]*X4[1]
J12 = dN1[1]*X1[2] + dN2[1]*X2[2] + dN3[1]*X3[2] + dN4[1]*X4[2]
# ... (all 9 components)
# Inverse Jacobian (manual 3x3 inverse)
detJ = J11*(J22*J33 - J23*J32) - J12*(J21*J33 - J23*J31) + ...
invJ11 = (J22*J33 - J23*J32) / detJ
# ... (all 9 components)
# Physical derivatives (manual matrix-vector multiply)
dN1_dx = (invJ11*dN1[1] + invJ12*dN1[2] + invJ13*dN1[3],
invJ21*dN1[1] + invJ22*dN1[2] + invJ23*dN1[3],
invJ31*dN1[1] + invJ32*dN1[2] + invJ33*dN1[3])
# ... (for all 4 nodes)
# Strain (manual tensor product, sum, symmetrize)
gradu11 = dN1_dx[1]*u1[1] + dN2_dx[1]*u2[1] + dN3_dx[1]*u3[1] + dN4_dx[1]*u4[1]
# ... (all 9 components)
ε = (gradu11,
gradu22,
gradu33,
(gradu12 + gradu21) / 2,
(gradu13 + gradu31) / 2,
(gradu23 + gradu32) / 2)
# Plasticity
state_old = states[elem_idx, 1] # 1 GP for Tet4
(σ, state_new, plastic) = return_mapping_gpu(ε, state_old, E, ν, σ_y)
states[elem_idx, 1] = state_new
# Forces (manual contraction)
f1 = (dN1_dx[1]*σ[1] + dN1_dx[2]*σ[4] + dN1_dx[3]*σ[5],
dN1_dx[1]*σ[4] + dN1_dx[2]*σ[2] + dN1_dx[3]*σ[6],
dN1_dx[1]*σ[5] + dN1_dx[2]*σ[6] + dN1_dx[3]*σ[3])
# ... (for all 4 nodes)
# Gauss weight * detJ
wdetJ = (1.0/6.0) * detJ
# Scatter to global (ATOMIC!)
CUDA.@atomic r_global[3*n1-2] += f1[1] * wdetJ
CUDA.@atomic r_global[3*n1-1] += f1[2] * wdetJ
CUDA.@atomic r_global[3*n1] += f1[3] * wdetJ
# ... (for all 4 nodes)
end
return nothing
end
```
## Next Steps
1. **Implement `element_residual_kernel!`** - The core GPU kernel (all manual indexing)
2. **Test kernel correctness** - Compare GPU vs CPU residual
3. **Wrap in matrix-free operator** - For GMRES
4. **Test GMRES with CuArrays** - Krylov.jl should handle it
5. **Complete Newton loop** - Anderson on CPU is fine
The key insight: **Everything except Anderson can stay on GPU!**
-77
View File
@@ -1,77 +0,0 @@
# Running GPU and MPI Demonstrations
## Prerequisites
These demonstrations require MPI and CUDA to be installed globally:
```bash
# MPI is already installed on your system
# CUDA is already installed on your system
```
The packages are loaded dynamically, so they don't need to be in `Project.toml`.
## Running the Demonstrations
### MPI Communication Test
Test data transfer between 2 MPI processes:
```bash
cd /home/juajukka/dev/JuliaFEM.jl
mpirun -np 2 julia benchmarks/gpu_mpi_demo.jl
```
Expected output:
- Rank 0 sends typed arrays to Rank 1
- Shows bytes transferred
- Validates data integrity
### GPU + MPI Combined Test
If CUDA GPU is available, runs full workflow:
```bash
mpirun -np 2 julia benchmarks/gpu_mpi_demo.jl
```
Expected output:
- Detects GPU (if available)
- Compiles type-stable kernel for GPU
- Executes on real hardware
- Transfers results via MPI
### Single-Process GPU Test
To test GPU without MPI:
```bash
julia benchmarks/gpu_only_demo.jl
```
## What Gets Demonstrated
1. **Type Stability Requirement**
- `Matrix{Float64}` transfers to GPU ✅
- `Dict{String,Any}` would FAIL GPU compilation ❌
2. **MPI Fast Transfer**
- Typed arrays: Fast buffer transfer (memcpy)
- Mixed types: Slow serialization (~100× slower)
3. **Real Hardware Execution**
- Actual CUDA kernel compilation and execution
- Actual MPI inter-process communication
- No mocks, no simulation
## Interpreting Results
Success indicators:
- ✅ "MPI communication successful" - Type-stable data transferred
- ✅ "GPU execution successful" - Kernel compiled and ran on GPU
- ✅ "Combined GPU+MPI workflow successful" - End-to-end validated
The key insight: The same type-stable patterns that give 9-92× CPU speedup also enable GPU execution and efficient MPI communication.
-212
View File
@@ -1,212 +0,0 @@
---
title: "GPU Assembly Proof-of-Concept"
date: 2025-11-10
status: "Working POC"
last_updated: 2025-11-10
tags: ["gpu", "proof-of-concept", "matrix-free", "elasticity"]
---
## Summary
**✅ PROOF OF CONCEPT COMPLETE!**
We have a working GPU assembly implementation that proves the entire finite element solve can stay on GPU with no escapes until the final result.
## What Works
### File: `demos/gpu_assembly_poc.jl`
**Features:**
- Element-parallel GPU kernel for 2D linear elasticity
- Quad4 elements with 2×2 Gauss quadrature
- Matrix-free Jacobian-vector product (finite difference on GPU)
- Complete Newton-Krylov loop on GPU
- Boundary condition enforcement
**Architecture:**
```
u0 (CPU) → GPU
[GPU Newton Loop]
- compute_residual_gpu!() # Element-parallel kernel
- compute_Jv_gpu!() # Finite difference Jv
- gmres() # Krylov.jl solver
- u .+= du # Update on GPU
u_final (CPU) ← GPU
```
**Validation:**
- ✅ GPU assembly matches CPU (relative error < 1e-15)
- ✅ Entire solve stays on GPU
- ✅ Only transfers: mesh data (once), u0 (in), u_final (out)
## Test Case
```julia
# Mesh: 10×10 Quad4 elements
Elements: 100
Nodes: 121
DOFs: 242
# Material: Steel
E = 200 GPa
ν = 0.3
# BC: Fixed left edge, displacement on right edge
```
## Performance
**Current Status (100 elements):**
- GPU assembly correct ✅
- Newton convergence: Issue (should converge in 1-2 iters for linear elasticity)
- Residual norm reduces but doesn't reach tolerance
**Known Issues:**
1. Finite difference epsilon may need tuning
2. GMRES tolerance may be too loose
3. BC enforcement could be improved
**Next:** Test larger problems (1K, 5K, 10K DOFs) to see GPU speedup
## Key Code Components
### GPU Kernel (Element-Parallel)
```julia
@cuda threads=256 blocks=n_blocks function elasticity_residual_kernel!(
r_global, u_global, elem_nodes, node_coords, E, ν
)
elem_id = threadIdx().x + (blockIdx().x - 1) * blockDim().x
# Compute element residual: r_elem = ∫ Bᵀ σ dV
for ip in 1:4
# Gauss quadrature, shape functions, strain, stress
...
end
# Atomic scatter to global residual
CUDA.@atomic r_global[dof] += r_elem[i]
end
```
### Matrix-Free Jv on GPU
```julia
function compute_Jv_gpu!(Jv, u, v, r0, ...)
u_perturbed = u .+ ε .* v # GPU vector operation
r_perturbed = compute_residual_gpu!(u_perturbed)
Jv .= (r_perturbed .- r0) ./ ε
end
```
## What This Proves
1.**GPU kernel correctness** - Assembly matches CPU to machine precision
2.**Matrix-free on GPU** - Jv computed via finite difference, all on GPU
3.**Resident data** - u, r stay on GPU entire solve
4.**Krylov.jl integration** - Works with GPU vectors
5.**Minimal transfers** - Only initial/final data moves
## Architecture Decisions Validated
From design documents (`docs/design/gpu_*.md`):
-**Element-parallel** works (atomic scatter acceptable for now)
-**Flat arrays** (CuMatrix, CuVector) - correct data structure
-**StaticArrays** for local operations (efficient)
-**No warp optimization needed yet** (proves concept first)
## Next Steps
### Short Term (Optimization)
1. **Debug Newton convergence** - Should be 1-2 iterations for linear elasticity
2. **Benchmark performance** - Test 1K, 5K, 10K DOFs
3. **Measure speedup** - Compare vs CPU assembly
4. **Profile kernels** - Memory bandwidth, atomic contention
### Medium Term (Integration)
1. **Create AssemblyState struct** - Proper data management
2. **GPU data structures** - to_gpu/to_cpu conversions
3. **Refactor into src/gpu/** - Proper module structure
4. **Multiple element types** - Hex8, Tet4, etc.
### Long Term (Research)
1. **Warp reduction** - Optimize atomics (32× reduction)
2. **Node-parallel kernel** - Your research idea
3. **Plasticity on GPU** - Material state updates
4. **Contact mechanics** - Node-based constraints
## Dependencies (Global)
```julia
using CUDA # GPU programming
using StaticArrays # Fast local arrays
using Krylov # Matrix-free solvers
```
**Note:** Installed globally, not in project environment per user request.
## Running the Demo
```bash
julia demos/gpu_assembly_poc.jl
```
**Output:**
```
======================================================================
GPU Assembly Proof-of-Concept
======================================================================
📐 Mesh:
Elements: 100 (Quad4)
Nodes: 121
DOFs: 242
🧪 Validating GPU vs CPU assembly...
Max absolute error: 6.4e-10
Relative error: 4.1e-16
✅ GPU assembly matches CPU!
🚀 Starting GPU Newton-Krylov solve...
Newton iter 1: ||r|| = 1.18e9
...
✅ PROOF OF CONCEPT COMPLETE!
======================================================================
```
## Lessons Learned
1. **CUDA.sync() → CUDA.synchronize()** - API naming
2. **Krylov.jl eltype warning** - Expected (operator wraps CuArrays)
3. **StaticArrays essential** - Fast local element operations
4. **Atomic scatter acceptable** - Not bottleneck yet for small problems
5. **BC enforcement** - Need to zero residual and du at fixed DOFs
## Comparison to Benchmark
From `benchmarks/matrix_free_gpu_benchmark.jl`:
- 3-4× GPU speedup at 5K-10K DOFs
- Matrix-free 3-8× faster than traditional Newton
**POC validates same architecture!**
## Design Documents
See comprehensive documentation:
1. `docs/design/gpu_assembly_architecture.md`
2. `docs/design/gpu_implementation_strategy.md`
3. `docs/design/gpu_kernel_comparison.md`
## Key Takeaway
> "Everything stays on GPU. No escapes until final result."
**Mission accomplished!**
The architecture is sound, the kernel is correct, and we've proven that GPU-resident matrix-free FEM is viable in Julia.
Now we optimize and integrate.
-307
View File
@@ -1,307 +0,0 @@
# Multi-GPU MPI Krylov Solver Demonstration
## Overview
This demonstration proves that type-stable nodal assembly enables distributed solving on multi-GPU systems using Krylov iterative methods. This is the complete workflow for modern scalable FEM.
## What This Demonstrates
### 1. Nodal Assembly Pattern
- **Row-by-row matrix construction**: Each rank assembles its local rows
- **`get_row()` abstraction**: Simulates nodal assembly from element contributions
- **Natural for contact mechanics**: Contact is inherently nodal
### 2. Distributed Computing
- **Domain decomposition**: 10×10 problem split across 2 MPI ranks
- **Each rank owns 5 nodes** (rows 1-5 and 6-10)
- **MPI collectives**: Global dot products via `Allreduce`, vector assembly via `Allgatherv`
### 3. Multi-GPU Execution
- **Local GPU per rank**: Each rank transfers data to its GPU
- **GPU matrix-vector products**: Computation on GPU, synchronization via MPI
- **Type-stable kernels**: GPU requires concrete types (no `Dict{String,Any}`)
### 4. Krylov Iterative Solver
- **Conjugate Gradient (CG)**: Matrix-free iterative solver
- **Distributed matvec**: Each rank computes `y_local = A_local * x_global`
- **Convergence**: 9 iterations to reach relative error < 1e-13
## Running the Demo
```bash
# With 2 MPI processes (recommended)
mpiexec -np 2 julia --project=. benchmarks/krylov_mpi_gpu_demo.jl
# With 4 processes (if you have 4 GPUs)
mpiexec -np 4 julia --project=. benchmarks/krylov_mpi_gpu_demo.jl
```
## Expected Output
```text
======================================================================
Multi-GPU MPI Krylov Solver Demonstration
======================================================================
Configuration:
MPI ranks: 2
CUDA available: true
Part 1: Generating Test Problem
----------------------------------------------------------------------
Problem size: 10×10 system
✓ Generated SPD matrix (condition number ≈ 3.45)
✓ Exact solution: x = [1, 2, 3, ..., 10]
Part 2: Nodal Assembly Pattern
----------------------------------------------------------------------
Rank 0: assembled 5 rows locally
✓ Nodal assembly complete (each rank has its partition)
Part 3: GPU Transfer
----------------------------------------------------------------------
Rank 0: transferred 440 bytes to GPU
✓ Each rank transferred local data to its GPU
Part 4: Distributed Matrix-Vector Product
----------------------------------------------------------------------
✓ Distributed matrix-vector product working
Part 5: Conjugate Gradient Solver
----------------------------------------------------------------------
Initial residual: 5.962915e+02
Iteration 1: residual = 7.590913e+01 (reduction: 87.27%)
Iteration 2: residual = 6.773957e+00 (reduction: 98.86%)
...
Iteration 9: residual = 1.835156e-11 (reduction: 100.00%)
✓ Converged in 9 iterations
Part 6: Verification
----------------------------------------------------------------------
Relative error: 7.731369e-14
✓ VERIFICATION PASSED
```
## Technical Details
### Problem Setup
- **Matrix**: 10×10 symmetric positive definite (SPD)
- **Condition number**: ~3.45 (well-conditioned)
- **Exact solution**: `x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]`
- **Right-hand side**: `b = A * x_exact`
### Partitioning
```text
Rank 0: owns nodes 1-5 (rows 1-5 of matrix)
Rank 1: owns nodes 6-10 (rows 6-10 of matrix)
```
### Distributed Matrix-Vector Product
Each rank:
1. Has local rows `A_local` (5×10 matrix)
2. Needs global vector `x_global` (10 entries)
3. Computes local result `y_local = A_local * x_global` (5 entries)
4. No communication needed during matvec (only during assembly of global vectors)
### CG Algorithm (Distributed)
```julia
# Initialization
r = b - A*x # Distributed: each rank has r_local
p = r # Search direction (needs to be global)
for iter = 1:maxiter
# Matrix-vector product (distributed)
Ap = A * p
# Global dot products (MPI Allreduce)
alpha = (r'*r) / (p'*Ap)
# Update solution and residual
x = x + alpha * p
r = r - alpha * Ap
# Check convergence
if ||r|| < tol
break
end
# Update search direction
beta = (r_new'*r_new) / (r_old'*r_old)
p = r + beta * p
end
```
### GPU Execution
- **CPU → GPU**: Transfer `A_local` (5×10 matrix) and `x_global` (10 vector)
- **GPU computation**: `y_local = A_local * x_global` (80 FLOPs)
- **GPU → CPU**: Transfer result `y_local` (5 entries)
- **Key requirement**: Type-stable data (`Matrix{Float64}`, not `Dict`)
### MPI Communication
- **`MPI.Allreduce`**: Sum scalar values across ranks (dot products)
- **`MPI.Allgatherv`**: Gather variable-length vectors from all ranks
- **Frequency**: Once per CG iteration (not during matvec)
## Performance Characteristics
### Communication Cost
- **Per CG iteration**:
- 2× `Allreduce` (scalar): ~O(log n_ranks) latency
- 2× `Allgatherv` (vector): ~O(N) bandwidth
- Total: Dominated by vector transfers, not latency
### Computation Cost
- **Per CG iteration**:
- 1× matvec: O(N²/n_ranks) FLOPs per rank
- 2× dot products: O(N/n_ranks) FLOPs per rank
- Total: O(N²/n_ranks) FLOPs
### Scaling
- **Weak scaling**: Problem size N increases with n_ranks → constant time
- **Strong scaling**: Fixed N, increase n_ranks → speedup until communication dominates
- **This demo**: Strong scaling (fixed N=10, tiny problem)
## Relevance to JuliaFEM
### Why This Matters
1. **Nodal assembly is natural for contact mechanics**
- Contact constraints are nodal (not element-based)
- Row-by-row assembly aligns with contact detection
- Streaming assembly: process nodes as they're detected
2. **Type stability is not optional**
- GPU kernels REQUIRE concrete types
- MPI benefits from typed buffers (no serialization)
- Dict-based fields CANNOT work in this workflow
3. **Matrix-free is the future**
- Don't assemble global matrix (memory scales as N²)
- Only need matvec operator (memory scales as N)
- Krylov methods only need matrix-vector products
4. **Distributed solving is achievable**
- Even small problems (10×10) work correctly
- Algorithm scales to millions of DOFs
- Same code works on 1 core, 1 GPU, or 100 GPUs
### Demonstrated Path: v0.5 → v1.0
**v0.5.1 (2019):**
- Element assembly → global matrix → direct solver
- Dict-based fields (type-unstable)
- Single-threaded CPU only
- Max ~100K DOF (memory limited)
**v1.0 (target):**
- Nodal assembly → matrix-free matvec → Krylov solver
- Type-stable fields (requirement)
- Multi-GPU + MPI (demonstrated here)
- Max ~10M DOF (computation limited)
### Next Steps
1. **Scale up problem size**: Test with N=1000, N=10000
2. **Add preconditioning**: Jacobi, ILU, AMG
3. **Real FEM integration**: Replace `get_row()` with actual assembly
4. **Performance profiling**: Measure communication vs computation ratio
5. **Strong scaling study**: Fix N, vary n_ranks, measure speedup
## Key Insights
### Type Stability Enables Everything
| Feature | Requires Type Stability? | Why? |
|---------|-------------------------|------|
| Fast CPU code | ✅ Yes | Avoid dispatch overhead (9-92× measured) |
| GPU execution | ✅ REQUIRED | Cannot compile kernels with abstract types |
| MPI transfer | ✅ Yes | Typed buffers avoid serialization |
| Krylov solvers | ✅ Yes | Matrix-free operators need concrete types |
| Threading | ✅ Yes | Race-free parallel access needs known layouts |
### Nodal Assembly Advantages
1. **Contact mechanics alignment**: Natural for nodal constraints
2. **Streaming assembly**: Process nodes incrementally (memory efficient)
3. **Domain decomposition**: Each rank owns nodes (clean partitioning)
4. **Matrix-free**: Never form global matrix (scalability)
### Krylov vs Direct Solvers
| Method | Memory | Time | Scalability | Notes |
|--------|--------|------|-------------|-------|
| Direct (LU) | O(N²) | O(N³) | Poor | v0.5.1 used this |
| Krylov (CG) | O(N) | O(N·iter) | Excellent | This demo |
| Krylov + Precond | O(N) | O(N·iter/√κ) | Best | Future work |
*κ = condition number, iter = iterations to converge*
## Validation
### Test Problem
- **Type**: Linear system `Ax = b`
- **Matrix**: 10×10 SPD, condition number ~3.45
- **Exact solution**: `x = [1, 2, 3, ..., 10]`
### Results
- **Converged in**: 9 iterations
- **Final residual**: 1.84 × 10⁻¹¹
- **Relative error**: 7.73 × 10⁻¹⁴
- **Status**: ✅ PASSED (error < 1e-6)
### Hardware
- **GPU**: NVIDIA RTX A2000 12GB (per rank)
- **MPI ranks**: 2
- **Data transferred**: 440 bytes per rank to GPU
- **GPU execution**: Verified working on both ranks
## References
### Conjugate Gradient Method
- Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain"
- Saad, "Iterative Methods for Sparse Linear Systems"
### Domain Decomposition
- Smith, Bjørstad, Gropp, "Domain Decomposition: Parallel Multilevel Methods for Elliptic PDEs"
### GPU Computing
- CUDA.jl documentation: https://cuda.juliagpu.org/
- Sanders, Kandrot, "CUDA by Example"
### MPI
- MPI.jl documentation: https://juliaparallel.org/MPI.jl/
- Gropp, Lusk, Skjellum, "Using MPI"
## Conclusion
This demonstration proves that:
1.**Type-stable nodal assembly works** on real hardware
2.**Multi-GPU execution is possible** with type-stable data
3.**MPI communication is efficient** with typed buffers
4.**Krylov solvers converge correctly** in distributed setting
5.**Solution accuracy is excellent** (relative error ~1e-13)
**The path forward is clear**: Type stability is the foundation, nodal assembly is the pattern, Krylov+MPI is the solver. This is how JuliaFEM v1.0 will scale to millions of DOFs.
-264
View File
@@ -1,264 +0,0 @@
---
title: "GPU POC: Tensors.jl Integration"
date: 2025-11-10
status: "Corrected Architecture"
last_updated: 2025-11-10
tags: ["gpu", "tensors", "architecture", "material-modeling"]
---
## The Problem
The initial GPU proof-of-concept (`gpu_assembly_poc.jl`) **ignored** the material modeling architecture established in `docs/book/material_modeling.md`.
**What was wrong:**
```julia
# ❌ OLD: Manual Voigt-like indexing
ε = SA[εxx, εyy, γxy] # Just a vector!
σ = C * ε # Matrix multiplication
r_elem[1] += (dN_dx[1] * σ[1] + dN_dy[1] * σ[3]) * factor # Manual indexing
```
**Problems:**
- No `SymmetricTensor` - just plain vectors
- No material API - hardcoded constitutive matrix
- Manual index arithmetic for stress components
- Doesn't match the established architecture!
## The Solution
**Corrected version** (`gpu_assembly_poc_tensors.jl`) uses proper Tensors.jl:
```julia
# ✅ NEW: Proper tensor operations
ε = SymmetricTensor{2,2}((εxx, γxy/2, εyy)) # Symmetric tensor!
σ = compute_stress_2d(material, ε) # Material API!
r_contrib = compute_B_transpose_sigma(dN_dx, dN_dy, σ) # Clean operations
```
**Advantages:**
-`SymmetricTensor{2,2}` for strain and stress (2D)
- ✅ Material API: `compute_stress(material, ε)`
- ✅ Follows `material_modeling.md` architecture
- ✅ GPU compatible (Tensors.jl works on CUDA!)
- ✅ Mathematics looks like equations
## Key Changes
### 1. Material Model Struct
```julia
struct LinearElastic
E::Float64
ν::Float64
end
@inline λ(mat::LinearElastic) = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
@inline μ(mat::LinearElastic) = mat.E / (2(1 + mat.ν))
```
### 2. Material API
```julia
@inline function compute_stress_2d(
material::LinearElastic,
ε::SymmetricTensor{2,2,T}
) where T
λ_val = T(λ(material))
μ_val = T(μ(material))
I = one(ε)
# Hooke's law: σ = λ·tr(ε)·I + 2μ·ε
σ = λ_val * tr(ε) * I + 2μ_val * ε
return σ
end
```
### 3. Strain Computation
```julia
@inline function compute_B_matrix_strain(dN_dx, dN_dy, u_elem)
"""Returns SymmetricTensor{2,2} for 2D strain"""
εxx = dN_dx[1] * u_elem[1] + ...
εyy = dN_dy[1] * u_elem[2] + ...
γxy = dN_dy[1] * u_elem[1] + dN_dx[1] * u_elem[2] + ...
# SymmetricTensor{2,2}: (ε11, ε12, ε22)
# Note: ε12 = γxy/2 (tensorial, not engineering shear)
return SymmetricTensor{2,2}((εxx, γxy/2, εyy))
end
```
### 4. Stress-to-Force Conversion
```julia
@inline function compute_B_transpose_sigma(dN_dx, dN_dy, σ::SymmetricTensor{2,2})
"""Compute Bᵀ·σ for element residual"""
# Extract stress components (automatic with Tensors.jl)
σxx = σ[1,1]
σyy = σ[2,2]
σxy = σ[1,2] # Symmetric, not engineering
# Nodal forces
r_elem = SA[
dN_dx[1] * σxx + dN_dy[1] * σxy, # Node 1, x
dN_dy[1] * σyy + dN_dx[1] * σxy, # Node 1, y
...
]
return r_elem
end
```
### 5. GPU Kernel
```julia
function elasticity_residual_kernel_tensors!(
r_global, u_global, elem_nodes, node_coords, E, ν
)
# Material model
material = LinearElastic(E, ν)
for ip in 1:4
# ...compute dN_dx, dN_dy...
# ✅ Tensor strain
ε = compute_B_matrix_strain(dN_dx, dN_dy, u_elem)
# ✅ Material API
σ = compute_stress_2d(material, ε)
# ✅ Clean force computation
r_contrib = compute_B_transpose_sigma(dN_dx, dN_dy, σ)
r_elem .+= r_contrib .* (w * det_J)
end
# Atomic scatter (same as before)
end
```
## Benefits
### 1. Extensibility
Adding new materials is **trivial**:
```julia
struct NeoHookean
C10::Float64
D1::Float64
end
@inline function compute_stress_2d(
material::NeoHookean,
ε::SymmetricTensor{2,2,T}
) where T
# Neo-Hookean stress computation
# Just define this function - kernel stays unchanged!
...
end
```
**GPU kernel doesn't change at all!** Dispatch handles it.
### 2. Plasticity Ready
```julia
struct VonMisesPlasticity
E::Float64
ν::Float64
σ_y::Float64 # Yield stress
end
struct PlasticState{T}
ε_p::SymmetricTensor{2,2,T} # Plastic strain
α::T # Hardening parameter
end
@inline function compute_stress_2d(
material::VonMisesPlasticity,
ε::SymmetricTensor{2,2,T},
state_old::PlasticState{T}
) where T
# Trial stress
ε_e = ε - state_old.ε_p
σ_trial = compute_stress_2d(LinearElastic(material.E, material.ν), ε_e)
# Check yield
σ_dev = dev(σ_trial) # Tensors.jl provides this!
σ_eq = (3/2 * σ_dev σ_dev) # von Mises stress
if σ_eq < material.σ_y
return σ_trial, state_old # Elastic
else
# Return mapping (closed-form for perfect plasticity)
...
end
end
```
**This is the architecture from `material_modeling.md`!**
### 3. Code Clarity
Compare old vs new for von Mises calculation:
```julia
# ❌ OLD (Voigt notation):
σ_dev = σ_vec - sum(σ_vec[1:3])/3 * [1,1,1,0,0,0]
σ_eq = (σ_dev[1]^2 + σ_dev[2]^2 + σ_dev[3]^2 +
2*(σ_dev[4]^2 + σ_dev[5]^2 + σ_dev[6]^2))
# ✅ NEW (Tensors.jl):
σ_dev = dev(σ)
σ_eq = (3/2 * σ_dev σ_dev)
```
**Mathematics looks like equations!**
## Current Status
### ✅ Working
- `gpu_assembly_poc_tensors.jl` runs on GPU
- Uses proper `SymmetricTensor{2,2}` for strain/stress
- Material API: `compute_stress_2d(material, ε)`
- Follows `material_modeling.md` architecture
- Extensible to new materials via dispatch
### ⚠️ Same Convergence Issue
Both versions have the same Newton convergence problem (doesn't converge for linear elasticity). This is a separate issue with:
- Finite difference epsilon size
- Boundary condition enforcement
- GMRES tolerance
**The kernel is correct** (same residual as CPU), convergence is secondary optimization.
## Files
- **Old (wrong):** `demos/gpu_assembly_poc.jl` - Manual indexing, no material API
- **New (correct):** `demos/gpu_assembly_poc_tensors.jl` - Proper Tensors.jl
- **Reference:** `docs/book/material_modeling.md` - Established architecture
## Next Steps
1.**Use Tensors.jl** - Done!
2. **Fix convergence** - Debug Newton/GMRES
3. **Add plasticity** - Implement `VonMisesPlasticity` material
4. **Benchmark** - Test 1K, 5K, 10K DOFs
5. **Integrate** - Move to `src/gpu/` proper architecture
## Key Takeaway
> "Always follow the established architecture in `docs/book/material_modeling.md`!"
The POC proved the GPU concept works, but **the second version proves it works with the correct architecture**.
---
**Lesson learned:** When user says "you forgot Tensors.jl and material_modeling.md", they're right! Always check design documents before coding.
-131
View File
@@ -1,131 +0,0 @@
# Assembly Strategy Comparison - Simple Example
#
# Demonstrates the modern Physics API for solving elasticity problems.
# Uses the CPU backend with element assembly.
using JuliaFEM
using LinearAlgebra
using Printf
println("="^70)
println("Assembly Comparison - Modern Physics API")
println("="^70)
# ============================================================================
# 1. Create Simple Mesh
# ============================================================================
println("\n[1] Creating mesh...")
# Simple 2-element beam (Hex8 elements)
nodes = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [2.0, 0.0, 0.0],
4 => [0.0, 1.0, 0.0],
5 => [1.0, 1.0, 0.0],
6 => [2.0, 1.0, 0.0],
7 => [0.0, 0.0, 1.0],
8 => [1.0, 0.0, 1.0],
9 => [2.0, 0.0, 1.0],
10 => [0.0, 1.0, 1.0],
11 => [1.0, 1.0, 1.0],
12 => [2.0, 1.0, 1.0]
)
connectivity_hex = [
(1, 2, 5, 4, 7, 8, 11, 10),
(2, 3, 6, 5, 8, 9, 12, 11)
]
n_nodes = length(nodes)
n_elements = length(connectivity_hex)
n_dofs = 3 * n_nodes
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs: $n_dofs")
# ============================================================================
# 2. Create Physics Problem
# ============================================================================
println("\n[2] Creating physics problem...")
physics = Physics(Elasticity, "simple beam", 3)
physics.properties.formulation = :continuum
physics.properties.finite_strain = false
# Create elements with new immutable API
elements = Element[]
for conn in connectivity_hex
# Extract node coordinates
X = [nodes[i] for i in conn]
# Create immutable element with all fields
element = Element(Hex8, conn,
fields=(geometry=X,
youngs_modulus=210e9, # Steel
poissons_ratio=0.3))
push!(elements, element)
end
add_elements!(physics, elements)
println(" Elements added: $(length(physics.body_elements))")
# ============================================================================
# 3. Apply Boundary Conditions
# ============================================================================
println("\n[3] Applying boundary conditions...")
# Fix left end (nodes 1, 4, 7, 10)
fixed_nodes = [1, 4, 7, 10]
add_dirichlet!(physics, fixed_nodes, [1, 2, 3], 0.0)
println(" Fixed nodes: $(length(fixed_nodes)) (all DOFs)")
println(" Total Dirichlet BCs: $(length(physics.bc_dirichlet.node_ids))")
# Note: External forces would be applied via Neumann BC or body forces
# For this simple demo, we solve with zero external loading
# ============================================================================
# 4. Solve with CPU Backend
# ============================================================================
println("\n[4] Solving with CPU backend...")
t_solve = @elapsed begin
sol = solve!(physics; backend=CPU(), tol=1e-6, max_iter=1000)
end
println(" Solve time: $(round(t_solve * 1000, digits=2)) ms")
println(" CG iterations: $(sol.cg_iterations)")
println(" Newton iterations: $(sol.newton_iterations)")
println(" Residual: $(sol.residual)")
println(" Max displacement: $(maximum(abs.(sol.u)) * 1000) mm")
# ============================================================================
# 5. Summary
# ============================================================================
println("\n" * "="^70)
println("SUMMARY")
println("="^70)
println("\nProblem:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs: $n_dofs")
println(" Fixed DOFs: $(3 * length(fixed_nodes))")
println("\nSolution:")
println(" Backend: CPU (element assembly + CG)")
println(" Solve time: $(round(t_solve * 1000, digits=2)) ms")
println(" CG iterations: $(sol.cg_iterations)")
println(" Newton iterations: $(sol.newton_iterations)")
println(" Residual: $(sol.residual)")
println(" Max displacement: $(maximum(abs.(sol.u)) * 1000) mm")
println("\n" * "="^70)
println("✓ Modern Physics API working on CPU!")
println("="^70)
-145
View File
@@ -1,145 +0,0 @@
"""
Demo: Solve cantilever beam problem with GPU
This demonstrates the complete workflow:
1. Generate mesh with Gmsh
2. Read mesh
3. Define material and boundary conditions
4. Solve on GPU
5. Visualize results
"""
# Load the GPU elasticity module
include(joinpath(@__DIR__, "..", "src", "gpu_elasticity.jl"))
using .GPUElasticity
using .GPUElasticity.GmshReader: get_surface_nodes
using Printf
function main()
println("\n" * "="^70)
println("Demo: Cantilever Beam on GPU")
println("="^70)
# Step 1: Generate mesh (if not exists)
mesh_file = joinpath(@__DIR__, "..", "test", "testdata", "cantilever_beam.msh")
if !isfile(mesh_file)
println("\nGenerating mesh...")
mkpath(dirname(mesh_file))
include(joinpath(@__DIR__, "..", "scripts", "generate_cantilever_mesh.jl"))
Base.invokelatest(generate_cantilever_mesh,
length=10.0,
width=1.0,
height=1.0,
mesh_size=0.5,
output_file=mesh_file
)
else
println("\nUsing existing mesh: $mesh_file")
end
# Step 2: Read mesh
println("\nReading mesh...")
mesh = read_gmsh_mesh(mesh_file)
# Step 3: Define boundary conditions
println("\nDefining boundary conditions...")
# Fixed end (X = 0)
fixed_nodes = get_surface_nodes(mesh, "FixedEnd")
println(" Fixed nodes: $(length(fixed_nodes))")
# Pressure surface (Z = max)
pressure_nodes = get_surface_nodes(mesh, "PressureSurface")
println(" Pressure nodes: $(length(pressure_nodes))")
# Step 4: Define material (steel)
material = ElasticMaterial(
210e9, # E = 210 GPa
0.3 # ν = 0.3
)
println("\nMaterial: Steel")
println(" E = $(material.E / 1e9) GPa")
println(" ν = $(material.ν)")
# Step 5: Define load (1 MPa pressure on top)
pressure = 1e6 # Pa
println("\nLoad: Pressure on top surface")
println(" Magnitude: $(pressure / 1e6) MPa")
# Step 6: Create physics
physics = ElasticityPhysics(
mesh,
material,
fixed_nodes,
pressure_nodes,
pressure
)
# Step 7: Solve on GPU
println("\n" * "="^70)
println("Solving on GPU...")
println("="^70)
u = solve_elasticity_gpu(physics, tol=1e-6, max_iter=1000)
# Step 8: Post-process results
println("\n" * "="^70)
println("Results:")
println("="^70)
n_nodes = size(mesh.nodes, 2)
# Displacement magnitudes
disp_mag = zeros(n_nodes)
for i in 1:n_nodes
ux = u[3*i-2]
uy = u[3*i-1]
uz = u[3*i]
disp_mag[i] = sqrt(ux^2 + uy^2 + uz^2)
end
println("\nDisplacement statistics:")
@printf(" Max: %.6e m\n", maximum(disp_mag))
@printf(" Min: %.6e m\n", minimum(disp_mag))
@printf(" Avg: %.6e m\n", sum(disp_mag) / n_nodes)
# Find node with max displacement
max_node = argmax(disp_mag)
x_max = mesh.nodes[1, max_node]
y_max = mesh.nodes[2, max_node]
z_max = mesh.nodes[3, max_node]
println("\nMax displacement location:")
@printf(" Node: %d\n", max_node)
@printf(" Position: (%.3f, %.3f, %.3f)\n", x_max, y_max, z_max)
@printf(" Displacement: (%.6e, %.6e, %.6e) m\n",
u[3*max_node-2], u[3*max_node-1], u[3*max_node])
# Analytical comparison
L = 10.0
width = 1.0
height = 1.0
I = width * height^3 / 12
q = pressure * width
w_analytical = q * L^4 / (8 * material.E * I)
println("\nComparison with beam theory:")
@printf(" Analytical: %.6e m\n", w_analytical)
@printf(" FEM: %.6e m\n", maximum(disp_mag))
@printf(" Error: %.2f%%\n",
abs(w_analytical - maximum(disp_mag)) / w_analytical * 100)
println("\n" * "="^70)
println("✅ Demo complete!")
println("="^70)
println("\nNext steps:")
println(" - Visualize with ParaView (export to VTK)")
println(" - Try different mesh sizes")
println(" - Add more complex loading")
println(" - Test with nonlinear materials")
println("="^70 * "\n")
end
main()
-183
View File
@@ -1,183 +0,0 @@
# Cantilever Beam - Simple Example Using Real JuliaFEM API
#
# Demonstrates the modern Physics API for solving elasticity problems.
# Shows both direct and iterative solvers on a realistic cantilever beam.
using JuliaFEM
using LinearAlgebra
using Printf
println("="^70)
println("Cantilever Beam - Simple Example")
println("="^70)
# ============================================================================
# 1. Create Mesh with Gmsh
# ============================================================================
println("\n[1] Generating mesh...")
using Gmsh: gmsh
gmsh.initialize()
gmsh.model.add("cantilever")
# Beam geometry: 10m × 1m × 1m
L, W, H = 10.0, 1.0, 1.0
lc = 1.5 # Mesh size
box = gmsh.model.occ.addBox(0, 0, 0, L, W, H)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.mesh.generate(3)
# Extract nodes and connectivity
node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
nodes = reshape(node_coords, 3, length(node_tags))
elem_types, _, elem_node_tags_vec = gmsh.model.mesh.getElements(3)
tet4_idx = findfirst(t -> t == 4, elem_types)
elem_node_tags = elem_node_tags_vec[tet4_idx]
connectivity = reshape(Int.(elem_node_tags), 4, :)
n_nodes = size(nodes, 2)
n_elements = size(connectivity, 2)
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs: $(3 * n_nodes)")
gmsh.finalize()
# ============================================================================
# 2. Create Physics Problem
# ============================================================================
println("\n[2] Setting up physics...")
# Create elasticity problem
physics = Physics(Elasticity, "cantilever", 3)
# Create elements with geometry and material properties
elements = Element[]
for e in 1:n_elements
conn = Tuple(connectivity[:, e])
element = Element(Tet4, conn)
# Set geometry
X = Dict(i => nodes[:, connectivity[i, e]] for i in 1:4)
update!(element, "geometry", X)
# Set material properties (steel)
update!(element, "youngs modulus", 210e9) # Pa
update!(element, "poissons ratio", 0.3)
push!(elements, element)
end
add_elements!(physics, elements)
# Boundary conditions
# Fixed: nodes at X=0
fixed_nodes = findall(x -> abs(x) < 1e-10, nodes[1, :])
add_dirichlet!(physics, fixed_nodes, [1, 2, 3], 0.0)
# Loaded: nodes at X=L (apply point loads)
loaded_nodes = findall(x -> abs(x - L) < 1e-10, nodes[1, :])
F_total = -1000.0 # Total force in Z direction
f_per_node = F_total / length(loaded_nodes)
println(" Fixed nodes: $(length(fixed_nodes))")
println(" Loaded nodes: $(length(loaded_nodes))")
println(" Force per node: $(f_per_node) N")
# Apply loads as body forces (workaround until Neumann BC works)
for element in elements
conn = get_connectivity(element)
for node_id in conn
if node_id in loaded_nodes
# This is a simplified approach - proper implementation would use Neumann BC
# For now, we'll assemble and apply loads manually
end
end
end
# ============================================================================
# 3. Solve with Direct Solver
# ============================================================================
println("\n[3] Assembling system...")
# Assemble using existing JuliaFEM infrastructure
problem = Problem(Elasticity, "body", 3)
add_elements!(problem, elements)
# Apply Dirichlet BCs
bc = Problem(Dirichlet, "fixed", 3, "displacement")
bc_elements = Element[]
for node in fixed_nodes
# Create point element for BC
bc_el = Element(Poi1, [node])
update!(bc_el, "geometry", Dict(node => nodes[:, node]))
update!(bc_el, "displacement 1", 0.0)
update!(bc_el, "displacement 2", 0.0)
update!(bc_el, "displacement 3", 0.0)
push!(bc_elements, bc_el)
end
add_elements!(bc, bc_elements)
# Apply loads
load = Problem(Elasticity, "load", 3)
load_elements = Element[]
for node in loaded_nodes
el = elements[findfirst(e -> node in get_connectivity(e), elements)]
# Add to existing element
update!(el, "displacement load 3", f_per_node)
end
# Assemble
t_assembly = @elapsed begin
assemble!(problem, 0.0)
assemble!(bc, 0.0)
end
println(" Assembly time: $(round(t_assembly, digits=4)) s")
# Solve
println("\n[4] Solving with direct solver...")
t_solve = @elapsed begin
K = problem.assembly.K
f = problem.assembly.f
# Apply BCs
eliminate_boundary_conditions!(problem, bc)
# Solve
K_full = Matrix(K)
f_full = Vector(f)
u = K_full \ f_full
end
println(" Solve time: $(round(t_solve, digits=4)) s")
println(" Total time: $(round(t_assembly + t_solve, digits=4)) s")
println(" Max displacement: $(maximum(abs.(u)) * 1000) mm")
# ============================================================================
# 5. Summary
# ============================================================================
println("\n" * "="^70)
println("SOLUTION SUMMARY")
println("="^70)
println("\nProblem size:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs: $(3 * n_nodes)")
println(" Fixed DOFs: $(3 * length(fixed_nodes))")
println("\nResults:")
println(" Max displacement: $(maximum(abs.(u)) * 1000) mm")
println(" Assembly time: $(round(t_assembly, digits=4)) s")
println(" Solve time: $(round(t_solve, digits=4)) s")
println(" Total time: $(round(t_assembly + t_solve, digits=4)) s")
println("\n" * "="^70)
-498
View File
@@ -1,498 +0,0 @@
# Cantilever Beam - Assembly Strategy Benchmark (RESEARCH/DEVELOPMENT)
#
# ⚠️ NOTE: This is a LOW-LEVEL benchmark for algorithm research!
# ⚠️ For USER-FACING examples, see:
# ⚠️ - demos/assembly_comparison_simple.jl (uses real Problem API)
# ⚠️ - demos/cantilever_gmsh_gpu.jl (uses Physics API + GPU)
# ⚠️ - examples/linear_static.jl (complete workflow)
#
# This file compares three assembly/solver combinations at the structure level:
# 1. Element assembly + Direct solver (baseline)
# 2. Element assembly + Iterative CG
# 3. Nodal assembly + Matrix-free CG (research)
#
# Uses element_assembly_structures.jl and nodal_assembly_structures.jl directly.
# Not intended as example of user-facing API!
using LinearAlgebra
using SparseArrays
using Tensors
using Printf
# Import our assembly structures
include("../src/element_assembly_structures.jl")
include("../src/nodal_assembly_structures.jl")
println("="^70)
println("Cantilever Beam - CPU Assembly Comparison")
println("="^70)
# ============================================================================
# 1. Generate Mesh with Gmsh
# ============================================================================
println("\n[1] Generating mesh with Gmsh...")
# Simple beam: L=10, W=1, H=1
# Target ~20 Tet4 elements
using Gmsh: gmsh
gmsh.initialize()
gmsh.model.add("cantilever")
# Geometry
lc = 1.5 # Characteristic length (controls mesh density)
L, W, H = 10.0, 1.0, 1.0
# Create box
box = gmsh.model.occ.addBox(0, 0, 0, L, W, H)
gmsh.model.occ.synchronize()
# Set mesh size
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
# Generate 3D mesh
gmsh.model.mesh.generate(3)
# Extract nodes
node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
n_nodes = length(node_tags)
nodes = reshape(node_coords, 3, n_nodes)
println(" Nodes: $n_nodes")
# Extract Tet4 elements (type 4)
elem_types, elem_tags_vec, elem_node_tags_vec = gmsh.model.mesh.getElements(3)
tet4_idx = findfirst(t -> t == 4, elem_types) # Type 4 = Tet4
if tet4_idx === nothing
error("No Tet4 elements found!")
end
elem_node_tags = elem_node_tags_vec[tet4_idx]
n_elements = div(length(elem_node_tags), 4)
connectivity = reshape(Int.(elem_node_tags), 4, n_elements)
println(" Elements: $n_elements")
println(" DOFs: $(3 * n_nodes)")
gmsh.finalize()
# ============================================================================
# 2. Material Properties and BCs
# ============================================================================
println("\n[2] Setting up problem...")
# Material (steel)
E = 210e9 # Young's modulus [Pa]
ν = 0.3 # Poisson's ratio
# Boundary conditions
# Fixed: nodes at X=0
fixed_nodes = findall(x -> abs(x) < 1e-10, nodes[1, :])
println(" Fixed nodes: $(length(fixed_nodes))")
# Loaded: nodes at X=L (free end)
loaded_nodes = findall(x -> abs(x - L) < 1e-10, nodes[1, :])
println(" Loaded nodes: $(length(loaded_nodes))")
# Applied force (total 1000 N downward, distributed)
F_total = -1000.0 # Negative Z direction
f_per_node = F_total / length(loaded_nodes)
println(" Force per node: $(f_per_node) N")
# ============================================================================
# 3. Compute Element Stiffness Matrices (Shared by all methods)
# ============================================================================
#
# NOTE: This is a simplified reference implementation for benchmarking.
# For production use, see src/problems_elasticity.jl which includes:
# - Geometric nonlinearity, finite strain
# - Plasticity and advanced material models
# - Surface tractions, body forces
# - Integration with Problem/Element API
#
# This demo focuses on assembly strategy comparison, not material complexity.
# ============================================================================
println("\n[3] Computing element stiffness matrices...")
# Elasticity tensor (isotropic)
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
function compute_tet4_stiffness(X::Matrix{Float64}, E::Float64, ν::Float64)
# X: 3×4 matrix of node coordinates
# Returns: 12×12 element stiffness matrix
# Shape function derivatives in parent element (constant for Tet4)
dN_dξ = [-1.0 -1.0 -1.0;
1.0 0.0 0.0;
0.0 1.0 0.0;
0.0 0.0 1.0]
# Jacobian: J = dX/dξ
J = X * dN_dξ # 3×3
detJ = det(J)
if detJ <= 0
error("Negative Jacobian determinant!")
end
# Shape function derivatives in physical space
dN_dx = dN_dξ / J # 4×3
# B matrix (strain-displacement): 6×12
B = zeros(6, 12)
for i in 1:4
B[1, 3i-2] = dN_dx[i, 1] # ∂u/∂x
B[2, 3i-1] = dN_dx[i, 2] # ∂v/∂y
B[3, 3i] = dN_dx[i, 3] # ∂w/∂z
B[4, 3i-2] = dN_dx[i, 2] # ∂u/∂y
B[4, 3i-1] = dN_dx[i, 1] # ∂v/∂x
B[5, 3i-1] = dN_dx[i, 3] # ∂v/∂z
B[5, 3i] = dN_dx[i, 2] # ∂w/∂y
B[6, 3i-2] = dN_dx[i, 3] # ∂u/∂z
B[6, 3i] = dN_dx[i, 1] # ∂w/∂x
end
# Elasticity matrix (Voigt notation)
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
D = [λ+2μ λ λ 0 0 0;
λ λ+2μ λ 0 0 0;
λ λ λ+2μ 0 0 0;
0 0 0 μ 0 0;
0 0 0 0 μ 0;
0 0 0 0 0 μ]
# Element stiffness: K_e = ∫ B^T D B dV = B^T D B * V
# For Tet4: V = detJ / 6
V = abs(detJ) / 6.0
K_e = (B' * D * B) * V
return K_e
end
# Compute all element matrices
K_elements = Vector{Matrix{Float64}}(undef, n_elements)
for e in 1:n_elements
conn = connectivity[:, e]
X_elem = nodes[:, conn]
K_elements[e] = compute_tet4_stiffness(X_elem, E, ν)
end
println(" Element stiffness matrices computed")
# ============================================================================
# 4. METHOD 1: Element Assembly + Direct Solver
# ============================================================================
println("\n" * "="^70)
println("METHOD 1: Element Assembly + Direct Solver (LU)")
println("="^70)
t1 = time()
# Assemble global system
n_dofs = 3 * n_nodes
assembly = ElementAssemblyData(n_dofs, Float64)
for e in 1:n_elements
conn = Tuple(connectivity[:, e])
gdofs = get_dof_indices(conn, 3)
contrib = ElementContribution(e, gdofs, K_elements[e],
zeros(12), zeros(12))
scatter_to_global!(assembly, contrib)
end
# Apply loads
for node in loaded_nodes
dof_z = 3 * node # Z component
assembly.f_ext_global[dof_z] = f_per_node
end
# Compute residual
compute_residual!(assembly)
# Apply Dirichlet BCs
fixed_dofs = Int[]
for node in fixed_nodes
append!(fixed_dofs, [3 * node - 2, 3 * node - 1, 3 * node])
end
apply_dirichlet_bc!(assembly, fixed_dofs, zeros(length(fixed_dofs)))
t_assembly_1 = time() - t1
println("Assembly time: $(round(t_assembly_1, digits=4)) s")
# Solve with direct solver (K_global is already CSC)
t_solve_1_start = time()
u1 = assembly.K_global \ assembly.r_global
t_solve_1 = time() - t_solve_1_start
# Compute final residual
r1 = matrix_vector_product(assembly, u1) - assembly.f_ext_global
r1_norm = norm(r1)
t_total_1 = time() - t1
println("Solve time: $(round(t_solve_1, digits=4)) s")
println("Total time: $(round(t_total_1, digits=4)) s")
println("Residual norm: $(r1_norm)")
println("Max displacement: $(maximum(abs.(u1)) * 1000) mm")
# ============================================================================
# 5. METHOD 2: Element Assembly + Iterative Solver (CG)
# ============================================================================
println("\n" * "="^70)
println("METHOD 2: Element Assembly + Iterative Solver (CG)")
println("="^70)
t2 = time()
# Reuse assembly from Method 1
t_assembly_2 = t_assembly_1 # Same assembly
# Conjugate Gradient solver
function cg_solve(A::ElementAssemblyData, b::Vector{Float64};
tol=1e-8, max_iter=1000)
n = length(b)
x = zeros(n)
r = b - matrix_vector_product(A, x)
p = copy(r)
rsold = dot(r, r)
for iter in 1:max_iter
Ap = matrix_vector_product(A, p)
α = rsold / dot(p, Ap)
x .+= α .* p
r .-= α .* Ap
rsnew = dot(r, r)
if sqrt(rsnew) < tol
return x, iter, sqrt(rsnew)
end
β = rsnew / rsold
p .= r .+ β .* p
rsold = rsnew
end
return x, max_iter, sqrt(rsold)
end
t_solve_2_start = time()
u2, cg_iters_2, cg_res_2 = cg_solve(assembly, assembly.r_global, tol=1e-8)
t_solve_2 = time() - t_solve_2_start
# Compute final residual
r2 = matrix_vector_product(assembly, u2) - assembly.f_ext_global
r2_norm = norm(r2)
t_total_2 = time() - t2
println("Assembly time: $(round(t_assembly_2, digits=4)) s")
println("Solve time: $(round(t_solve_2, digits=4)) s")
println("Total time: $(round(t_total_2, digits=4)) s")
println("CG iterations: $cg_iters_2")
println("CG residual: $(cg_res_2)")
println("Residual norm: $(r2_norm)")
println("Max displacement: $(maximum(abs.(u2)) * 1000) mm")
println("Difference from Method 1: $(norm(u1 - u2))")
# ============================================================================
# 6. METHOD 3: Nodal Assembly + Iterative Solver (Matrix-Free CG)
# ============================================================================
println("\n" * "="^70)
println("METHOD 3: Nodal Assembly + Matrix-Free Iterative Solver")
println("="^70)
t3 = time()
# Build node-to-elements map (convert matrix to vector of tuples)
conn_tuples = [Tuple(connectivity[:, e]) for e in 1:n_elements]
node_map = NodeToElementsMap(conn_tuples)
# For each node, precompute 3×3 stiffness blocks with all coupling nodes
# This is the "spider" pattern
struct NodalAssemblyData
node_map::NodeToElementsMap
K_elements::Vector{Matrix{Float64}}
connectivity::Matrix{Int}
conn_tuples::Vector{NTuple{4,Int}} # Store tuple version too
n_nodes::Int
n_dofs::Int
end
nodal_data = NodalAssemblyData(node_map, K_elements, connectivity,
conn_tuples, n_nodes, n_dofs)
# Matrix-vector product using nodal assembly
function nodal_matvec!(w::Vector{Float64}, v::Vector{Float64},
data::NodalAssemblyData)
fill!(w, 0.0)
for node_i in 1:data.n_nodes
# Get spider nodes (all nodes coupled to node_i)
spider = get_node_spider(data.node_map, node_i, data.conn_tuples)
w_local = zeros(3)
for node_j in spider
# Sum contributions from all elements containing both nodes
K_block_ij = zeros(3, 3)
for elem_info in data.node_map.node_to_elements[node_i]
elem_idx = elem_info.element_id
conn = data.conn_tuples[elem_idx]
# Check if node_j is in this element
local_j = findfirst(==(node_j), conn)
if local_j !== nothing
local_i = findfirst(==(node_i), conn)
K_e = data.K_elements[elem_idx]
# Extract 3×3 block
for α in 1:3, β in 1:3
K_block_ij[α, β] += K_e[3*(local_i-1)+α, 3*(local_j-1)+β]
end
end
end
# Apply to displacement
v_j = v[3*(node_j-1)+1:3*node_j]
w_local .+= K_block_ij * v_j
end
# Write to global
w[3*(node_i-1)+1:3*node_i] .= w_local
end
return w
end
t_assembly_3 = time() - t3
println("Nodal map construction: $(round(t_assembly_3, digits=4)) s")
# Build RHS (same as before)
f_ext = zeros(n_dofs)
for node in loaded_nodes
dof_z = 3 * node
f_ext[dof_z] = f_per_node
end
# CG with matrix-free matvec
function cg_solve_nodal(data::NodalAssemblyData, b::Vector{Float64},
fixed_dofs::Vector{Int};
tol=1e-8, max_iter=1000)
n = length(b)
x = zeros(n)
# Apply BC to initial guess
x[fixed_dofs] .= 0.0
# Compute initial residual
Ax = zeros(n)
nodal_matvec!(Ax, x, data)
Ax[fixed_dofs] .= 0.0 # Zero out fixed DOFs
r = b - Ax
r[fixed_dofs] .= 0.0
p = copy(r)
rsold = dot(r, r)
for iter in 1:max_iter
Ap = zeros(n)
nodal_matvec!(Ap, p, data)
Ap[fixed_dofs] .= 0.0
α = rsold / dot(p, Ap)
x .+= α .* p
r .-= α .* Ap
rsnew = dot(r, r)
if sqrt(rsnew) < tol
return x, iter, sqrt(rsnew)
end
β = rsnew / rsold
p .= r .+ β .* p
rsold = rsnew
end
return x, max_iter, sqrt(rsold)
end
t_solve_3_start = time()
u3, cg_iters_3, cg_res_3 = cg_solve_nodal(nodal_data, f_ext, fixed_dofs,
tol=1e-8)
t_solve_3 = time() - t_solve_3_start
# Compute final residual
w3 = zeros(n_dofs)
nodal_matvec!(w3, u3, nodal_data)
r3 = w3 - f_ext
r3_norm = norm(r3)
t_total_3 = time() - t3
println("Solve time: $(round(t_solve_3, digits=4)) s")
println("Total time: $(round(t_total_3, digits=4)) s")
println("CG iterations: $cg_iters_3")
println("CG residual: $(cg_res_3)")
println("Residual norm: $(r3_norm)")
println("Max displacement: $(maximum(abs.(u3)) * 1000) mm")
println("Difference from Method 1: $(norm(u1 - u3))")
# ============================================================================
# 7. Summary Comparison
# ============================================================================
println("\n" * "="^70)
println("SUMMARY COMPARISON")
println("="^70)
println("\nProblem Size:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs: $n_dofs")
println(" Fixed DOFs: $(length(fixed_dofs))")
println(" Free DOFs: $(n_dofs - length(fixed_dofs))")
println("\n" * "-"^70)
println(@sprintf("%-40s %10s %10s %10s", "Method", "Assembly", "Solve", "Total"))
println("-"^70)
println(@sprintf("%-40s %9.4fs %9.4fs %9.4fs",
"1. Element + Direct (LU)", t_assembly_1, t_solve_1, t_total_1))
println(@sprintf("%-40s %9.4fs %9.4fs %9.4fs",
"2. Element + Iterative (CG, $cg_iters_2 iter)",
t_assembly_2, t_solve_2, t_total_2))
println(@sprintf("%-40s %9.4fs %9.4fs %9.4fs",
"3. Nodal + Iterative (CG, $cg_iters_3 iter)",
t_assembly_3, t_solve_3, t_total_3))
println("-"^70)
println("\nAccuracy (vs Method 1):")
println(@sprintf(" Method 2 error: %.3e", norm(u1 - u2)))
println(@sprintf(" Method 3 error: %.3e", norm(u1 - u3)))
println("\nSpeedup vs Method 1:")
println(@sprintf(" Method 2: %.2fx", t_total_1 / t_total_2))
println(@sprintf(" Method 3: %.2fx", t_total_1 / t_total_3))
println("\nMax Displacement:")
println(@sprintf(" Method 1: %.6f mm", maximum(abs.(u1)) * 1000))
println(@sprintf(" Method 2: %.6f mm", maximum(abs.(u2)) * 1000))
println(@sprintf(" Method 3: %.6f mm", maximum(abs.(u3)) * 1000))
println("\n" * "="^70)
println("All methods complete!")
println("="^70)
-383
View File
@@ -1,383 +0,0 @@
# Cantilever Beam Demo - GPU Physics with Gmsh Mesh
#
# Uses Gmsh to generate a proper mesh with ~20 quadratic tetrahedrons (Tet10)
# Demonstrates:
# - Gmsh mesh generation (programmatic)
# - Quadratic elements (Tet10, Tri6)
# - Immutable element API
# - GPU solver with Physics{Elasticity}
using JuliaFEM
using CUDA
using Tensors
println("="^60)
println("Cantilever Beam - GPU Physics with Gmsh Mesh")
println("="^60)
# Include new GPU physics module
include("../src/gpu_physics_elasticity.jl")
using .GPUElasticityPhysics
# ============================================================================
# 1. Generate Mesh with Gmsh
# ============================================================================
println("\n[1] Generating mesh with Gmsh...")
# Check if gmsh is available
if !success(`which gmsh`)
error("Gmsh not found! Install with: sudo apt install gmsh (Linux) or brew install gmsh (Mac)")
end
# Create Gmsh script
geo_file = "/tmp/cantilever_beam.geo"
msh_file = "/tmp/cantilever_beam.msh"
open(geo_file, "w") do f
write(
f,
"""
// Cantilever beam geometry
// Dimensions: L=4.0, W=1.0, H=1.0
SetFactory("OpenCASCADE");
// Create box
Box(1) = {0, 0, 0, 4.0, 1.0, 1.0};
// Define physical groups
Physical Volume("body") = {1};
// Fixed end (X=0)
Physical Surface("fixed") = {1};
// Free end (X=4) - not needed for load
Physical Surface("free") = {2};
// Top surface (Z=1) for pressure load
Physical Surface("pressure") = {6};
// Mesh settings
Mesh.CharacteristicLengthMin = 0.3;
Mesh.CharacteristicLengthMax = 0.5;
Mesh.ElementOrder = 2; // Quadratic elements (Tet10)
Mesh.Algorithm3D = 4; // Frontal Delaunay
// Generate 3D mesh
Mesh 3;
"""
)
end
# Run Gmsh
println(" Running Gmsh...")
run(`gmsh $geo_file -3 -o $msh_file -format msh2`)
println(" ✓ Mesh generated: $msh_file")
# ============================================================================
# 2. Read Mesh (Simple Parser for MSH2 Format)
# ============================================================================
println("\n[2] Reading mesh...")
function read_msh2(filename::String)
nodes = Dict{Int,Vector{Float64}}()
elements = Dict{String,Vector{Vector{Int}}}()
physical_groups = Dict{String,Vector{Vector{Int}}}() # name => elements
open(filename) do f
section = ""
while !eof(f)
line = strip(readline(f))
if line == "\$Nodes"
section = "nodes"
n_nodes = parse(Int, readline(f))
for _ in 1:n_nodes
parts = split(readline(f))
node_id = parse(Int, parts[1])
x, y, z = parse(Float64, parts[2]), parse(Float64, parts[3]), parse(Float64, parts[4])
nodes[node_id] = [x, y, z]
end
elseif line == "\$Elements"
section = "elements"
n_elements = parse(Int, readline(f))
for _ in 1:n_elements
parts = split(readline(f))
elem_id = parse(Int, parts[1])
elem_type = parse(Int, parts[2])
n_tags = parse(Int, parts[3])
# Read tags
physical_tag = n_tags >= 1 ? parse(Int, parts[4]) : 0
# Connectivity starts after tags
offset = 4 + n_tags
# Element type: 4=Tet4, 11=Tet10, 2=Tri3, 9=Tri6
if elem_type == 11 # Tet10
connectivity = [parse(Int, parts[i]) for i in offset:offset+9]
if !haskey(elements, "Tet10")
elements["Tet10"] = Vector{Int}[]
end
push!(elements["Tet10"], connectivity)
# Store by physical group
group_name = "body_$physical_tag"
if !haskey(physical_groups, group_name)
physical_groups[group_name] = Vector{Int}[]
end
push!(physical_groups[group_name], connectivity)
elseif elem_type == 9 # Tri6
connectivity = [parse(Int, parts[i]) for i in offset:offset+5]
if !haskey(elements, "Tri6")
elements["Tri6"] = Vector{Int}[]
end
push!(elements["Tri6"], connectivity)
# Store by physical group (surface elements)
group_name = "surface_$physical_tag"
if !haskey(physical_groups, group_name)
physical_groups[group_name] = Vector{Int}[]
end
push!(physical_groups[group_name], connectivity)
elseif elem_type == 2 # Tri3
connectivity = [parse(Int, parts[i]) for i in offset:offset+2]
if !haskey(elements, "Tri3")
elements["Tri3"] = Vector{Int}[]
end
push!(elements["Tri3"], connectivity)
# Store by physical group
group_name = "surface_$physical_tag"
if !haskey(physical_groups, group_name)
physical_groups[group_name] = Vector{Int}[]
end
push!(physical_groups[group_name], connectivity)
end
end
elseif line == "\$EndNodes" || line == "\$EndElements"
section = ""
end
end
end
return nodes, elements, physical_groups
end
node_dict, elem_dict, physical_groups = read_msh2(msh_file)
println("\nPhysical groups found:")
for (name, elems) in physical_groups
println(" $name: $(length(elems)) elements")
end
# Convert nodes to matrix
n_nodes = length(node_dict)
node_ids = sort(collect(keys(node_dict)))
node_map = Dict(id => i for (i, id) in enumerate(node_ids))
nodes = zeros(3, n_nodes)
for (id, i) in node_map
nodes[:, i] = node_dict[id]
end
println(" Nodes: $n_nodes")
println(" Tet10 elements: $(length(get(elem_dict, "Tet10", [])))")
println(" Tri6 elements: $(length(get(elem_dict, "Tri6", [])))")
# ============================================================================
# 3. Create Body Elements (Tet10)
# ============================================================================
println("\n[3] Creating body elements...")
body_elements = Element[]
for conn_gmsh in elem_dict["Tet10"]
# Renumber nodes (Gmsh → 1-based sequential)
conn = [node_map[id] for id in conn_gmsh]
# Extract coordinates
X_elem = nodes[:, conn]
# Create immutable element
el = Element(Tet10, Lagrange{Tet10,2}, tuple(conn...);
fields=(geometry=X_elem,
youngs_modulus=210e9, # Steel
poissons_ratio=0.3))
push!(body_elements, el)
end
println(" Body elements created: $(length(body_elements))")
# ============================================================================
# 4. Create Physics{Elasticity}
# ============================================================================
println("\n[4] Creating Physics{Elasticity}...")
physics = Physics(Elasticity, "cantilever beam", 3)
physics.properties.formulation = :continuum
physics.properties.finite_strain = false
add_elements!(physics, body_elements)
println(" Physics: $(physics.name)")
println(" Elements in physics: $(length(physics.body_elements))")
# ============================================================================
# 5. Add Boundary Conditions
# ============================================================================
println("\n[5] Adding boundary conditions...")
# Dirichlet BC: Use "fixed" physical surface (tag 1)
# Extract all unique nodes from fixed surface elements
fixed_nodes_set = Set{Int}()
if haskey(physical_groups, "surface_1") # Physical tag 1 = "fixed"
for conn_gmsh in physical_groups["surface_1"]
for node_id in conn_gmsh
push!(fixed_nodes_set, node_map[node_id])
end
end
end
fixed_nodes = sort(collect(fixed_nodes_set))
println(" Fixing $(length(fixed_nodes)) nodes from 'fixed' surface (all DOFs)")
for node in fixed_nodes
add_dirichlet!(physics, [node], [1, 2, 3], 0.0)
end
# Neumann BC: Use "pressure" physical surface (tag 3)
println(" Applying pressure load on 'pressure' surface...")
pressure = -1e6 # -1 MPa in -Z direction
traction = Vec{3}((0.0, 0.0, pressure))
n_pressure_surfaces = 0
if haskey(physical_groups, "surface_3") # Physical tag 3 = "pressure"
for conn_gmsh in physical_groups["surface_3"]
conn = [node_map[id] for id in conn_gmsh]
X_surf = nodes[:, conn]
surf_el = Element(Tri6, Lagrange{Tri6,2}, tuple(conn...);
fields=(geometry=X_surf,))
add_neumann!(physics, surf_el, traction)
n_pressure_surfaces += 1
end
end
println(" Dirichlet BCs: $(length(physics.bc_dirichlet.node_ids)) nodes")
println(" Neumann BCs: $n_pressure_surfaces surface elements")
# ============================================================================
# 6. Solve on GPU
# ============================================================================
println("\n[6] Solving on GPU...")
println(" DOFs: $(3 * n_nodes)")
println(" Initializing GPU data...")
result = solve_elasticity_gpu!(
physics;
time=0.0,
tol=1e-6,
max_iter=2000 # More DOFs, might need more iterations
)
println("\n" * "="^60)
println("SOLUTION")
println("="^60)
println(" Iterations: $(result.iterations)")
println(" Residual: $(result.residual)")
# ============================================================================
# 7. Post-Process Results
# ============================================================================
println("\n[7] Post-processing...")
u = result.u
n_dofs = length(u)
# Extract displacement components
u_x = u[1:3:end]
u_y = u[2:3:end]
u_z = u[3:3:end]
# Compute magnitude
u_mag = sqrt.(u_x .^ 2 + u_y .^ 2 + u_z .^ 2)
println("\nDisplacement Statistics:")
println(" Max |u|: $(maximum(u_mag) * 1000) mm")
println(" Max u_x: $(maximum(abs.(u_x)) * 1000) mm")
println(" Max u_y: $(maximum(abs.(u_y)) * 1000) mm")
println(" Max u_z: $(maximum(abs.(u_z)) * 1000) mm")
# Find max displacement location
max_idx = argmax(u_mag)
println("\nMax displacement at node $max_idx:")
println(" Location: $(nodes[:, max_idx])")
println(" u = [$(u_x[max_idx]*1000), $(u_y[max_idx]*1000), $(u_z[max_idx]*1000)] mm")
# ============================================================================
# 8. Validate
# ============================================================================
println("\n[8] Validation...")
# Find free end nodes from "free" physical surface (tag 2)
free_end_nodes_set = Set{Int}()
if haskey(physical_groups, "surface_2") # Physical tag 2 = "free"
for conn_gmsh in physical_groups["surface_2"]
for node_id in conn_gmsh
push!(free_end_nodes_set, node_map[node_id])
end
end
end
free_end_nodes = collect(free_end_nodes_set)
free_end_disp = maximum(u_mag[free_end_nodes])
fixed_end_disp = maximum(u_mag[fixed_nodes])
println(" Free end max |u|: $(free_end_disp * 1000) mm")
println(" Fixed end max |u|: $(fixed_end_disp * 1000) mm")
if fixed_end_disp < 1e-10
println(" ✓ Fixed end has zero displacement (good!)")
else
println(" ✗ Fixed end displacement > 0 (bad!)")
end
if free_end_disp > 1e-6
println(" ✓ Free end has non-zero displacement (good!)")
else
println(" ✗ Free end displacement ≈ 0 (bad!)")
end
println("\n" * "="^60)
println("Demo complete!")
println("="^60)
# ============================================================================
# Summary
# ============================================================================
println("\nMesh Statistics:")
println(" Nodes: $n_nodes")
println(" DOFs: $(3 * n_nodes)")
println(" Tet10 elements: $(length(body_elements))")
println(" Tri6 surface elements: $n_pressure_surfaces")
println()
println("Performance:")
println(" CG iterations: $(result.iterations)")
println(" Final residual: $(result.residual)")
println(" GPU solver: Pure device code (no CPU-GPU transfer)")
-261
View File
@@ -1,261 +0,0 @@
# Cantilever Beam Demo - Pure GPU Physics{Elasticity}
#
# Demonstrates new API:
# - Physics{Elasticity} (not Problem)
# - Elements store geometry via update!(element, "geometry", nodes)
# - BCs via add_dirichlet! and add_neumann!
# - Pure GPU solve with solve_elasticity_gpu!
# IMPORTANT: Load CUDA before JuliaFEM for GPU backend support
using CUDA
using JuliaFEM
using Tensors
println("="^60)
println("Cantilever Beam - Backend-Transparent Demo")
println("="^60)
# ============================================================================
# 1. Create Mesh (Simple hand-coded for demo)
# ============================================================================
println("\n[1] Creating mesh...")
# 2×1×1 beam with 2 Tet4 elements (minimal example)
nodes = [
# Beam: L=2, W=1, H=1
0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 # X
0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 # Y
0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 # Z
]
# Tet4 connectivity (two elements spanning the beam)
tet_connectivity = [
[1, 2, 3, 5], # Element 1
[2, 3, 4, 6], # Element 2
[3, 4, 5, 7], # Element 3
[4, 5, 6, 8], # Element 4
]
println(" Nodes: $(size(nodes, 2))")
println(" Elements: $(length(tet_connectivity))")
# ============================================================================
# 2. Create Body Elements with Geometry and Material
# ============================================================================
println("\n[2] Creating body elements...")
body_elements = Element[]
for conn in tet_connectivity
# Create element with geometry and material (immutable API)
X_elem = nodes[:, conn]
# Legacy API: topology + connectivity (infers Lagrange{Tet4,1})
el = Element(Tet4, conn;
fields=(geometry=X_elem,
youngs_modulus=210e9, # 210 GPa (steel)
poissons_ratio=0.3))
push!(body_elements, el)
end
println(" Body elements created: $(length(body_elements))")
# ============================================================================
# 3. Create Physics{Elasticity}
# ============================================================================
println("\n[3] Creating Physics{Elasticity}...")
physics = Physics(Elasticity, "cantilever beam", 3)
physics.properties.formulation = :continuum
physics.properties.finite_strain = false
add_elements!(physics, body_elements)
println(" Physics: $(physics.name)")
println(" Elements in physics: $(length(physics.body_elements))")
# ============================================================================
# 4. Add Boundary Conditions
# ============================================================================
println("\n[4] Adding boundary conditions...")
# Dirichlet BC: Fix nodes at X=0 (left end)
fixed_nodes = [1, 3, 5, 7] # Nodes with X=0
println(" Fixing nodes: $fixed_nodes (all DOFs)")
for node in fixed_nodes
add_dirichlet!(physics, [node], [1, 2, 3], 0.0) # Fix u_x, u_y, u_z = 0
end
# Neumann BC: Pressure on top surface (Z=1)
# Surface: nodes [5, 6, 7, 8] form two triangles
println(" Applying pressure load on top surface...")
pressure = -1e6 # -1 MPa in -Z direction
traction = Vec{3}((0.0, 0.0, pressure))
# Top surface triangles
top_surface_tris = [
[5, 6, 7],
[6, 7, 8]
]
for tri_conn in top_surface_tris
X_surf = nodes[:, tri_conn]
# Legacy API: topology + connectivity (infers basis)
surf_el = Element(Tri3, tri_conn;
fields=(geometry=X_surf,))
add_neumann!(physics, surf_el, traction)
end
println(" Dirichlet BCs: $(length(physics.bc_dirichlet.node_ids)) nodes")
println(" Neumann BCs: $(length(physics.bc_neumann.surface_elements)) surfaces")
# ============================================================================
# 5. Solve (Backend Transparent!)
# ============================================================================
println("\n[5] Solving...")
println(" Using GPU backend (CPU backend not yet implemented)")
# Unified solve! - backend transparent!
# Explicitly request GPU since CPU backend is just a stub
result = solve!(
physics;
backend=GPU(), # Use GPU (CPU backend coming in Phase 2)
time=0.0,
tol=1e-6,
max_iter=1000
)
println("\n" * "="^60)
println("SOLUTION")
println("="^60)
if result.newton_iterations == 1
# Linear problem (converged in 1 Newton iteration)
println(" Problem type: LINEAR elasticity")
println(" Newton iterations: 1 (linear problem - no iterations needed)")
println(" CG iterations: $(result.cg_iterations)")
println(" Final residual: $(result.residual)")
println()
println(" Note: For LINEAR problems, Newton converges in 1 iteration.")
println(" CG iterations shown are from the single linear solve.")
else
# Nonlinear problem (multiple Newton iterations)
println(" Problem type: NONLINEAR elasticity")
println(" Solver: Inexact Newton-Krylov (SIMULTANEOUS solving!)")
println()
println(" Newton iterations: $(result.newton_iterations)")
println(" Total CG iterations: $(result.cg_iterations)")
println(" Avg CG per Newton: $(round(result.cg_iterations / result.newton_iterations, digits=1))")
println(" Final residual: $(result.residual)")
println()
println(" Iteration history (Newton + CG solved SIMULTANEOUSLY):")
for (i, (cg_i, R_i, η_i)) in enumerate(result.history)
println(" Newton $i: CG=$cg_i, ||R||=$(round(R_i, sigdigits=3)), η=$(round(η_i, digits=3))")
end
println()
println(" ✓ Linear systems solved INEXACTLY (adaptive tolerance)")
println(" ✓ Newton and Krylov iterations INTERLEAVED")
println(" ✓ Much faster than nested loops!")
end
println()
println(" Solve time: $(round(result.solve_time, digits=3)) seconds")
# ============================================================================
# 6. Post-Process Results
# ============================================================================
println("\n[6] Post-processing...")
u = result.u
n_nodes = div(length(u), 3)
# Extract displacement components
u_x = u[1:3:end]
u_y = u[2:3:end]
u_z = u[3:3:end]
# Compute magnitude
u_mag = sqrt.(u_x .^ 2 + u_y .^ 2 + u_z .^ 2)
println("\nDisplacement Statistics:")
println(" Max |u|: $(maximum(u_mag) * 1000) mm")
println(" Max u_x: $(maximum(abs.(u_x)) * 1000) mm")
println(" Max u_y: $(maximum(abs.(u_y)) * 1000) mm")
println(" Max u_z: $(maximum(abs.(u_z)) * 1000) mm")
# Find max displacement location
max_idx = argmax(u_mag)
println("\nMax displacement at node $max_idx:")
println(" Location: $(nodes[:, max_idx])")
println(" u = [$(u_x[max_idx]*1000), $(u_y[max_idx]*1000), $(u_z[max_idx]*1000)] mm")
# ============================================================================
# 7. Validate (Simple Check)
# ============================================================================
println("\n[7] Validation...")
# For cantilever with pressure load, expect:
# - Free end (X=2) has largest displacement
# - Fixed end (X=0) has zero displacement
# - Deflection primarily in -Z direction
free_end_nodes = [2, 4, 6, 8] # X=2
fixed_end_nodes = [1, 3, 5, 7] # X=0
free_end_disp = maximum(u_mag[free_end_nodes])
fixed_end_disp = maximum(u_mag[fixed_end_nodes])
println(" Free end max |u|: $(free_end_disp * 1000) mm")
println(" Fixed end max |u|: $(fixed_end_disp * 1000) mm")
if fixed_end_disp < 1e-10
println(" ✓ Fixed end has zero displacement (good!)")
else
println(" ✗ Fixed end displacement > 0 (bad!)")
end
if free_end_disp > 1e-6
println(" ✓ Free end has non-zero displacement (good!)")
else
println(" ✗ Free end displacement ≈ 0 (bad!)")
end
println("\n" * "="^60)
println("Demo complete!")
println("="^60)
# ============================================================================
# Summary
# ============================================================================
println("\nBackend-Transparent API Summary:")
println(" 1. Create elements: Element(Topology, connectivity; fields=(...))")
println(" 2. Create Physics: physics = Physics(Elasticity, name, dimension)")
println(" 3. Add elements: add_elements!(physics, body_elements)")
println(" 4. Add Dirichlet: add_dirichlet!(physics, node_ids, components, value)")
println(" 5. Add Neumann: add_neumann!(physics, surface_element, traction)")
println(" 6. Solve: result = solve!(physics) # ← Backend automatic!")
println()
println("Key Features:")
println(" ✓ Immutable elements with fields at construction")
println(" ✓ Backend transparency - same code for CPU or GPU")
println(" ✓ Automatic backend selection (GPU if CUDA available)")
println(" ✓ Matrix-free CG solver")
println(" ✓ No 'GPU' in user code!")
println()
println("Backend Selection:")
println(" solve!(physics) # Auto (GPU if available, else CPU)")
println(" solve!(physics; backend=GPU()) # Force GPU")
println(" solve!(physics; backend=CPU(8)) # Force CPU with 8 threads")
-606
View File
@@ -1,606 +0,0 @@
"""
GPU Assembly Proof-of-Concept
==============================
Minimal working example of matrix-free Newton-Krylov on GPU for 2D linear elasticity.
Goal: Prove that entire solve can stay on GPU with no escapes until final result.
Architecture:
- Element-parallel kernel (one thread per element)
- Simple atomic scatter (no warp optimization yet)
- Hardcoded Quad4 elements, 2x2 Gauss quadrature
- Matrix-free Jacobian-vector product via finite difference
- GMRES from Krylov.jl
- Everything stays on GPU during Newton loop
Status: PROOF OF CONCEPT - focus on correctness, optimize later
"""
using CUDA
using LinearAlgebra
using Krylov
# ============================================================================
# Mesh Generation: Simple rectangular mesh
# ============================================================================
function generate_rectangle_mesh(nx::Int, ny::Int, Lx::Float64, Ly::Float64)
"""Generate structured Quad4 mesh for rectangle [0,Lx] × [0,Ly]"""
# Node coordinates
n_nodes = (nx + 1) * (ny + 1)
coords = zeros(n_nodes, 2)
node_id = 1
for j in 0:ny
for i in 0:nx
coords[node_id, 1] = i * Lx / nx
coords[node_id, 2] = j * Ly / ny
node_id += 1
end
end
# Element connectivity (counterclockwise from lower-left)
n_elements = nx * ny
connectivity = zeros(Int32, n_elements, 4)
elem_id = 1
for j in 0:(ny-1)
for i in 0:(nx-1)
n1 = i + j * (nx + 1) + 1
n2 = (i + 1) + j * (nx + 1) + 1
n3 = (i + 1) + (j + 1) * (nx + 1) + 1
n4 = i + (j + 1) * (nx + 1) + 1
connectivity[elem_id, :] = [n1, n2, n3, n4]
elem_id += 1
end
end
return coords, connectivity
end
# ============================================================================
# GPU Kernel: Compute residual for linear elasticity
# ============================================================================
# Gauss quadrature points and weights (2x2 for Quad4)
const GAUSS_POINTS = SA[
SA[-0.5773502691896257, -0.5773502691896257],
SA[0.5773502691896257, -0.5773502691896257],
SA[0.5773502691896257, 0.5773502691896257],
SA[-0.5773502691896257, 0.5773502691896257]
]
const GAUSS_WEIGHTS = SA[1.0, 1.0, 1.0, 1.0]
@inline function shape_functions_quad4(ξ, η)
"""Quad4 shape functions at (ξ, η) ∈ [-1,1]²"""
return SA[
0.25*(1-ξ)*(1-η),
0.25*(1+ξ)*(1-η),
0.25*(1+ξ)*(1+η),
0.25*(1-ξ)*(1+η)
]
end
@inline function shape_derivatives_quad4(ξ, η)
"""Quad4 shape function derivatives: dN/dξ and dN/dη"""
dN_dξ = SA[
-0.25*(1-η),
0.25*(1-η),
0.25*(1+η),
-0.25*(1+η)
]
dN_dη = SA[
-0.25*(1-ξ),
-0.25*(1+ξ),
0.25*(1+ξ),
0.25*(1-ξ)
]
return dN_dξ, dN_dη
end
@inline function compute_jacobian_2d(dN_dξ, dN_dη, x_coords, y_coords)
"""Compute 2D Jacobian matrix: J = [dx/dξ dx/dη; dy/dξ dy/dη]"""
dx_dξ = sum(dN_dξ[i] * x_coords[i] for i in 1:4)
dx_dη = sum(dN_dη[i] * x_coords[i] for i in 1:4)
dy_dξ = sum(dN_dξ[i] * y_coords[i] for i in 1:4)
dy_dη = sum(dN_dη[i] * y_coords[i] for i in 1:4)
return SA[dx_dξ dy_dξ; dx_dη dy_dη] # Note: transposed for correct layout
end
@inline function constitutive_matrix_plane_strain(E, ν)
"""Plane strain constitutive matrix"""
factor = E / ((1 + ν) * (1 - 2ν))
return SA[
factor*(1-ν) factor*ν 0.0;
factor*ν factor*(1-ν) 0.0;
0.0 0.0 factor*(1-2ν)/2
]
end
# Main GPU kernel
function elasticity_residual_kernel!(
r_global::CuDeviceVector{T},
u_global::CuDeviceVector{T},
elem_nodes::CuDeviceMatrix{Int32},
node_coords::CuDeviceMatrix{T},
E::T,
ν::T
) where T
"""
Compute residual r = ∫ Bᵀ σ dV for linear elasticity.
One thread per element (element-parallel).
Uses atomic scatter for shared DOF contributions.
"""
elem_id = threadIdx().x + (blockIdx().x - 1) * blockDim().x
if elem_id > size(elem_nodes, 1)
return
end
# Get element nodes
n1, n2, n3, n4 = elem_nodes[elem_id, 1], elem_nodes[elem_id, 2],
elem_nodes[elem_id, 3], elem_nodes[elem_id, 4]
# Get node coordinates
x_coords = SA[node_coords[n1, 1], node_coords[n2, 1],
node_coords[n3, 1], node_coords[n4, 1]]
y_coords = SA[node_coords[n1, 2], node_coords[n2, 2],
node_coords[n3, 2], node_coords[n4, 2]]
# Get element DOFs (8 DOFs: 2 per node)
u_elem = SA[
u_global[2*n1-1], u_global[2*n1],
u_global[2*n2-1], u_global[2*n2],
u_global[2*n3-1], u_global[2*n3],
u_global[2*n4-1], u_global[2*n4]
]
# Constitutive matrix
C = constitutive_matrix_plane_strain(E, ν)
# Accumulate element residual
r_elem = MVector{8,T}(zeros(8))
# Loop over integration points
for ip in 1:4
ξ, η = GAUSS_POINTS[ip]
w = GAUSS_WEIGHTS[ip]
# Shape function derivatives
dN_dξ, dN_dη = shape_derivatives_quad4(ξ, η)
# Jacobian and its inverse
J = compute_jacobian_2d(dN_dξ, dN_dη, x_coords, y_coords)
det_J = J[1, 1] * J[2, 2] - J[1, 2] * J[2, 1]
inv_J = SA[J[2, 2] -J[1, 2]; -J[2, 1] J[1, 1]] / det_J
# Physical derivatives: [dN/dx; dN/dy] = inv(J) * [dN/dξ; dN/dη]
dN_dx = SA[
inv_J[1, 1]*dN_dξ[1]+inv_J[1, 2]*dN_dη[1],
inv_J[1, 1]*dN_dξ[2]+inv_J[1, 2]*dN_dη[2],
inv_J[1, 1]*dN_dξ[3]+inv_J[1, 2]*dN_dη[3],
inv_J[1, 1]*dN_dξ[4]+inv_J[1, 2]*dN_dη[4]
]
dN_dy = SA[
inv_J[2, 1]*dN_dξ[1]+inv_J[2, 2]*dN_dη[1],
inv_J[2, 1]*dN_dξ[2]+inv_J[2, 2]*dN_dη[2],
inv_J[2, 1]*dN_dξ[3]+inv_J[2, 2]*dN_dη[3],
inv_J[2, 1]*dN_dξ[4]+inv_J[2, 2]*dN_dη[4]
]
# B-matrix for strain-displacement (3×8)
# ε = [εxx, εyy, γxy]ᵀ = B * u_elem
# B = [dN1/dx 0 dN2/dx 0 dN3/dx 0 dN4/dx 0 ]
# [0 dN1/dy 0 dN2/dy 0 dN3/dy 0 dN4/dy]
# [dN1/dy dN1/dx dN2/dy dN2/dx dN3/dy dN3/dx dN4/dy dN4/dx]
# Compute strain: ε = B * u_elem
εxx = dN_dx[1] * u_elem[1] + dN_dx[2] * u_elem[3] +
dN_dx[3] * u_elem[5] + dN_dx[4] * u_elem[7]
εyy = dN_dy[1] * u_elem[2] + dN_dy[2] * u_elem[4] +
dN_dy[3] * u_elem[6] + dN_dy[4] * u_elem[8]
γxy = dN_dy[1] * u_elem[1] + dN_dx[1] * u_elem[2] +
dN_dy[2] * u_elem[3] + dN_dx[2] * u_elem[4] +
dN_dy[3] * u_elem[5] + dN_dx[3] * u_elem[6] +
dN_dy[4] * u_elem[7] + dN_dx[4] * u_elem[8]
ε = SA[εxx, εyy, γxy]
# Stress: σ = C * ε
σ = C * ε
# Add to element residual: r_elem += Bᵀ * σ * w * det(J)
factor = w * det_J
r_elem[1] += (dN_dx[1] * σ[1] + dN_dy[1] * σ[3]) * factor
r_elem[2] += (dN_dy[1] * σ[2] + dN_dx[1] * σ[3]) * factor
r_elem[3] += (dN_dx[2] * σ[1] + dN_dy[2] * σ[3]) * factor
r_elem[4] += (dN_dy[2] * σ[2] + dN_dx[2] * σ[3]) * factor
r_elem[5] += (dN_dx[3] * σ[1] + dN_dy[3] * σ[3]) * factor
r_elem[6] += (dN_dy[3] * σ[2] + dN_dx[3] * σ[3]) * factor
r_elem[7] += (dN_dx[4] * σ[1] + dN_dy[4] * σ[3]) * factor
r_elem[8] += (dN_dy[4] * σ[2] + dN_dx[4] * σ[3]) * factor
end
# Scatter to global residual (ATOMIC - multiple elements share nodes)
CUDA.@atomic r_global[2*n1-1] += r_elem[1]
CUDA.@atomic r_global[2*n1] += r_elem[2]
CUDA.@atomic r_global[2*n2-1] += r_elem[3]
CUDA.@atomic r_global[2*n2] += r_elem[4]
CUDA.@atomic r_global[2*n3-1] += r_elem[5]
CUDA.@atomic r_global[2*n3] += r_elem[6]
CUDA.@atomic r_global[2*n4-1] += r_elem[7]
CUDA.@atomic r_global[2*n4] += r_elem[8]
return nothing
end
# ============================================================================
# GPU Assembly Functions
# ============================================================================
function compute_residual_gpu!(
r_gpu::CuVector{T},
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T
) where T
"""Launch GPU kernel to compute residual"""
n_elements = size(elem_nodes_gpu, 1)
threads = 256
blocks = cld(n_elements, threads)
# Zero out residual
fill!(r_gpu, zero(T))
# Launch kernel
@cuda threads = threads blocks = blocks elasticity_residual_kernel!(
r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν
)
CUDA.synchronize()
return nothing
end
function compute_Jv_gpu!(
Jv_gpu::CuVector{T},
u_gpu::CuVector{T},
v_gpu::CuVector{T},
r0_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
ε::T=T(1e-7)
) where T
"""
Compute matrix-free Jacobian-vector product: Jv ≈ [R(u + εv) - R(u)] / ε
Everything stays on GPU!
"""
# Perturb u
u_perturbed = u_gpu .+ ε .* v_gpu # GPU vector operation
# Compute residual at perturbed state
r_perturbed = CUDA.zeros(T, length(u_gpu))
compute_residual_gpu!(r_perturbed, u_perturbed, elem_nodes_gpu, coords_gpu, E, ν)
# Finite difference approximation
Jv_gpu .= (r_perturbed .- r0_gpu) ./ ε
return nothing
end
# ============================================================================
# Matrix-Free Operator for Krylov.jl
# ============================================================================
struct GPUMatrixFreeOperator{T}
u::CuVector{T}
r0::CuVector{T}
elem_nodes::CuMatrix{Int32}
coords::CuMatrix{T}
E::T
ν::T
n::Int
end
function Base.size(op::GPUMatrixFreeOperator)
return (op.n, op.n)
end
function LinearAlgebra.mul!(Jv, op::GPUMatrixFreeOperator{T}, v) where T
"""Matrix-vector product for Krylov.jl"""
v_gpu = CuVector{T}(v)
Jv_gpu = CuVector{T}(undef, length(v))
compute_Jv_gpu!(Jv_gpu, op.u, v_gpu, op.r0, op.elem_nodes, op.coords, op.E, op.ν)
copyto!(Jv, Array(Jv_gpu))
return Jv
end
# ============================================================================
# GPU Newton-Krylov Solver
# ============================================================================
function solve_newton_krylov_gpu!(
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
fixed_dofs::Vector{Int};
max_iter::Int=20,
tol::T=T(1e-8),
gmres_tol::T=T(1e-6),
verbose::Bool=true
) where T
"""
Solve nonlinear elasticity problem using Newton-Krylov on GPU.
ENTIRE LOOP STAYS ON GPU - no escapes until convergence!
"""
n_dofs = length(u_gpu)
r_gpu = CUDA.zeros(T, n_dofs)
for iter in 1:max_iter
# Compute residual on GPU
compute_residual_gpu!(r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν)
# Enforce BC: zero out residual at fixed DOFs
r_cpu_temp = Array(r_gpu)
r_cpu_temp[fixed_dofs] .= 0.0
copyto!(r_gpu, r_cpu_temp)
# Check convergence (small data transfer for convergence check)
r_norm = CUDA.norm(r_gpu)
if verbose
println(" Newton iter $iter: ||r|| = $r_norm")
end
if r_norm < tol
if verbose
println(" ✅ Converged in $iter iterations")
end
return iter
end
# Matrix-free operator
op = GPUMatrixFreeOperator(u_gpu, r_gpu, elem_nodes_gpu, coords_gpu, E, ν, n_dofs)
# GMRES solve: J * du = -r
r_cpu = Array(-r_gpu)
du_cpu, stats = gmres(op, r_cpu, atol=gmres_tol, rtol=0.0, verbose=0)
if !stats.solved
@warn "GMRES did not converge at Newton iteration $iter"
end
# Enforce BC: zero out du at fixed DOFs
du_cpu[fixed_dofs] .= 0.0
# Update solution (transfer du back to GPU)
du_gpu = CuVector{T}(du_cpu)
u_gpu .+= du_gpu
end
@warn "Newton did not converge in $max_iter iterations"
return max_iter
end
# ============================================================================
# CPU Reference Implementation (for validation)
# ============================================================================
function compute_residual_cpu!(
r::Vector{T},
u::Vector{T},
elem_nodes::Matrix{Int32},
coords::Matrix{T},
E::T,
ν::T
) where T
"""CPU reference implementation"""
fill!(r, zero(T))
C = constitutive_matrix_plane_strain(E, ν)
for elem_id in 1:size(elem_nodes, 1)
n1, n2, n3, n4 = elem_nodes[elem_id, :]
x_coords = SA[coords[n1, 1], coords[n2, 1], coords[n3, 1], coords[n4, 1]]
y_coords = SA[coords[n1, 2], coords[n2, 2], coords[n3, 2], coords[n4, 2]]
u_elem = SA[
u[2*n1-1], u[2*n1],
u[2*n2-1], u[2*n2],
u[2*n3-1], u[2*n3],
u[2*n4-1], u[2*n4]
]
r_elem = MVector{8,T}(zeros(8))
for ip in 1:4
ξ, η = GAUSS_POINTS[ip]
w = GAUSS_WEIGHTS[ip]
dN_dξ, dN_dη = shape_derivatives_quad4(ξ, η)
J = compute_jacobian_2d(dN_dξ, dN_dη, x_coords, y_coords)
det_J = J[1, 1] * J[2, 2] - J[1, 2] * J[2, 1]
inv_J = SA[J[2, 2] -J[1, 2]; -J[2, 1] J[1, 1]] / det_J
dN_dx = SA[
inv_J[1, 1]*dN_dξ[1]+inv_J[1, 2]*dN_dη[1],
inv_J[1, 1]*dN_dξ[2]+inv_J[1, 2]*dN_dη[2],
inv_J[1, 1]*dN_dξ[3]+inv_J[1, 2]*dN_dη[3],
inv_J[1, 1]*dN_dξ[4]+inv_J[1, 2]*dN_dη[4]
]
dN_dy = SA[
inv_J[2, 1]*dN_dξ[1]+inv_J[2, 2]*dN_dη[1],
inv_J[2, 1]*dN_dξ[2]+inv_J[2, 2]*dN_dη[2],
inv_J[2, 1]*dN_dξ[3]+inv_J[2, 2]*dN_dη[3],
inv_J[2, 1]*dN_dξ[4]+inv_J[2, 2]*dN_dη[4]
]
εxx = dN_dx[1] * u_elem[1] + dN_dx[2] * u_elem[3] +
dN_dx[3] * u_elem[5] + dN_dx[4] * u_elem[7]
εyy = dN_dy[1] * u_elem[2] + dN_dy[2] * u_elem[4] +
dN_dy[3] * u_elem[6] + dN_dy[4] * u_elem[8]
γxy = dN_dy[1] * u_elem[1] + dN_dx[1] * u_elem[2] +
dN_dy[2] * u_elem[3] + dN_dx[2] * u_elem[4] +
dN_dy[3] * u_elem[5] + dN_dx[3] * u_elem[6] +
dN_dy[4] * u_elem[7] + dN_dx[4] * u_elem[8]
ε = SA[εxx, εyy, γxy]
σ = C * ε
factor = w * det_J
r_elem[1] += (dN_dx[1] * σ[1] + dN_dy[1] * σ[3]) * factor
r_elem[2] += (dN_dy[1] * σ[2] + dN_dx[1] * σ[3]) * factor
r_elem[3] += (dN_dx[2] * σ[1] + dN_dy[2] * σ[3]) * factor
r_elem[4] += (dN_dy[2] * σ[2] + dN_dx[2] * σ[3]) * factor
r_elem[5] += (dN_dx[3] * σ[1] + dN_dy[3] * σ[3]) * factor
r_elem[6] += (dN_dy[3] * σ[2] + dN_dx[3] * σ[3]) * factor
r_elem[7] += (dN_dx[4] * σ[1] + dN_dy[4] * σ[3]) * factor
r_elem[8] += (dN_dy[4] * σ[2] + dN_dx[4] * σ[3]) * factor
end
r[2*n1-1] += r_elem[1]
r[2*n1] += r_elem[2]
r[2*n2-1] += r_elem[3]
r[2*n2] += r_elem[4]
r[2*n3-1] += r_elem[5]
r[2*n3] += r_elem[6]
r[2*n4-1] += r_elem[7]
r[2*n4] += r_elem[8]
end
return nothing
end
# ============================================================================
# Main Demo
# ============================================================================
function main()
println("\n" * "="^70)
println("GPU Assembly Proof-of-Concept")
println("="^70)
# Problem setup
nx, ny = 10, 10 # 10×10 mesh = 100 elements, 121 nodes, 242 DOFs
Lx, Ly = 1.0, 1.0
E, ν = 200e9, 0.3 # Steel properties
println("\n📐 Mesh:")
println(" Elements: $(nx*ny) (Quad4)")
println(" Nodes: $((nx+1)*(ny+1))")
println(" DOFs: $(2*(nx+1)*(ny+1))")
# Generate mesh
coords, connectivity = generate_rectangle_mesh(nx, ny, Lx, Ly)
n_dofs = 2 * size(coords, 1)
println("\n🔧 Material:")
println(" Young's modulus: $(E/1e9) GPa")
println(" Poisson's ratio: $ν")
# Apply boundary conditions: fix left edge (x=0)
# and apply displacement on right edge (x=Lx)
fixed_dofs = Int[]
for node_id in 1:size(coords, 1)
if coords[node_id, 1] < 1e-10 # Left edge
push!(fixed_dofs, 2 * node_id - 1) # Fix x-displacement
push!(fixed_dofs, 2 * node_id) # Fix y-displacement
end
end
# Initial guess (small random perturbation)
u0 = randn(n_dofs) * 1e-6
# Apply Dirichlet BC: set fixed DOFs to zero
u0[fixed_dofs] .= 0.0
# Apply displacement BC on right edge (small tension)
for node_id in 1:size(coords, 1)
if abs(coords[node_id, 1] - Lx) < 1e-10 # Right edge
u0[2*node_id-1] = 0.001 # 1mm displacement in x
end
end
println("\n🔒 Boundary conditions:")
println(" Fixed DOFs: $(length(fixed_dofs))")
println(" Applied displacement: 1mm tension on right edge")
# Transfer to GPU
println("\n📤 Transferring data to GPU...")
elem_nodes_gpu = CuArray{Int32}(connectivity)
coords_gpu = CuArray{Float64}(coords)
u_gpu = CuArray{Float64}(u0)
println(" elem_nodes: $(size(elem_nodes_gpu))")
println(" coords: $(size(coords_gpu))")
println(" u: $(size(u_gpu))")
# Validate GPU assembly vs CPU
println("\n🧪 Validating GPU vs CPU assembly...")
r_cpu = zeros(n_dofs)
r_gpu = CUDA.zeros(Float64, n_dofs)
compute_residual_cpu!(r_cpu, u0, connectivity, coords, E, ν)
compute_residual_gpu!(r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν)
r_gpu_cpu = Array(r_gpu)
max_error = maximum(abs.(r_gpu_cpu .- r_cpu))
rel_error = max_error / (maximum(abs.(r_cpu)) + 1e-10)
println(" Max absolute error: $max_error")
println(" Relative error: $rel_error")
if rel_error < 1e-10
println(" ✅ GPU assembly matches CPU!")
else
println(" ❌ GPU assembly does NOT match CPU!")
return
end
# Solve using GPU Newton-Krylov
println("\n🚀 Starting GPU Newton-Krylov solve...")
println(" (Everything stays on GPU until convergence)")
u_gpu_solve = copy(u_gpu)
n_iter = solve_newton_krylov_gpu!(
u_gpu_solve, elem_nodes_gpu, coords_gpu, E, ν, fixed_dofs,
max_iter=20, tol=1e-8, gmres_tol=1e-6, verbose=true
)
# Transfer final solution back
u_final = Array(u_gpu_solve)
println("\n📊 Results:")
println(" Newton iterations: $n_iter")
println(" Final ||u||: $(norm(u_final))")
println(" Min displacement: $(minimum(u_final))")
println(" Max displacement: $(maximum(u_final))")
println("\n✅ PROOF OF CONCEPT COMPLETE!")
println(" - GPU assembly kernel works")
println(" - Matrix-free Jv on GPU works")
println(" - Newton loop stays on GPU (only u0 in, u_final out)")
println("="^70 * "\n")
end
# Run demo
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
-468
View File
@@ -1,468 +0,0 @@
"""
GPU Assembly Proof-of-Concept (Tensors.jl Version)
===================================================
Minimal working example of matrix-free Newton-Krylov on GPU for 2D linear elasticity.
**CORRECTED VERSION** using proper Tensors.jl architecture from material_modeling.md
Key changes from v1:
- Uses SymmetricTensor{2,2} for strain and stress (2D)
- Proper material API: compute_stress(material, ε, state, Δt)
- No Voigt notation, no manual indexing
- Mathematics looks like equations!
Goal: Prove that entire solve can stay on GPU with no escapes until final result.
"""
using CUDA
using LinearAlgebra
using Tensors # ✅ Using Tensors.jl for all tensor operations!
using Krylov
# ============================================================================
# Material Model (following material_modeling.md)
# ============================================================================
"""
Linear elastic (Hookean) material model for plane strain.
Stateless: σ depends only on current ε, no history.
"""
struct LinearElastic
E::Float64 # Young's modulus [Pa]
ν::Float64 # Poisson's ratio [-]
end
# Lamé parameters
@inline λ(mat::LinearElastic) = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
@inline μ(mat::LinearElastic) = mat.E / (2(1 + mat.ν))
"""
Compute stress for 2D plane strain using Tensors.jl.
Returns (σ, 𝔻, state_new) following unified material API.
"""
@inline function compute_stress_2d(
material::LinearElastic,
ε::SymmetricTensor{2,2,T}
) where T
λ_val = T(λ(material))
μ_val = T(μ(material))
# Identity tensor
I = one(ε)
# Hooke's law: σ = λ·tr(ε)·I + 2μ·ε
σ = λ_val * tr(ε) * I + 2μ_val * ε
return σ
end
# ============================================================================
# Mesh Generation
# ============================================================================
function generate_rectangle_mesh(nx::Int, ny::Int, Lx::Float64, Ly::Float64)
"""Generate structured Quad4 mesh for rectangle [0,Lx] × [0,Ly]"""
# Node coordinates
n_nodes = (nx + 1) * (ny + 1)
coords = zeros(n_nodes, 2)
node_id = 1
for j in 0:ny
for i in 0:nx
coords[node_id, 1] = i * Lx / nx
coords[node_id, 2] = j * Ly / ny
node_id += 1
end
end
# Element connectivity (counterclockwise from lower-left)
n_elements = nx * ny
connectivity = zeros(Int32, n_elements, 4)
elem_id = 1
for j in 0:(ny-1)
for i in 0:(nx-1)
n1 = i + j * (nx + 1) + 1
n2 = (i + 1) + j * (nx + 1) + 1
n3 = (i + 1) + (j + 1) * (nx + 1) + 1
n4 = i + (j + 1) * (nx + 1) + 1
connectivity[elem_id, :] = [n1, n2, n3, n4]
elem_id += 1
end
end
return coords, connectivity
end
# ============================================================================
# GPU Kernel using Tensors.jl
# ============================================================================
# Gauss quadrature (2x2 for Quad4)
const GAUSS_POINTS_2D = SA[
SA[-0.5773502691896257, -0.5773502691896257],
SA[0.5773502691896257, -0.5773502691896257],
SA[0.5773502691896257, 0.5773502691896257],
SA[-0.5773502691896257, 0.5773502691896257]
]
const GAUSS_WEIGHTS_2D = SA[1.0, 1.0, 1.0, 1.0]
@inline function shape_derivatives_quad4(ξ, η)
"""Quad4 shape function derivatives: dN/dξ and dN/dη"""
dN_dξ = SA[-0.25*(1-η), 0.25*(1-η), 0.25*(1+η), -0.25*(1+η)]
dN_dη = SA[-0.25*(1-ξ), -0.25*(1+ξ), 0.25*(1+ξ), 0.25*(1-ξ)]
return dN_dξ, dN_dη
end
@inline function compute_jacobian_quad4(dN_dξ, dN_dη, x_coords, y_coords)
"""Compute 2D Jacobian as Tensor{2,2}"""
dx_dξ = sum(dN_dξ[i] * x_coords[i] for i in 1:4)
dx_dη = sum(dN_dη[i] * x_coords[i] for i in 1:4)
dy_dξ = sum(dN_dξ[i] * y_coords[i] for i in 1:4)
dy_dη = sum(dN_dη[i] * y_coords[i] for i in 1:4)
# Return as Tensor (not SMatrix) for proper inv() and det()
return Tensor{2,2}((dx_dξ, dy_dξ, dx_dη, dy_dη))
end
@inline function compute_B_matrix_strain(dN_dx, dN_dy, u_elem)
"""
Compute strain from B-matrix and displacements using Tensors.jl.
Returns SymmetricTensor{2,2} for 2D strain.
"""
# ε = [εxx εxy] where εxy = (∂ux/∂y + ∂uy/∂x)/2
# [εxy εyy]
εxx = dN_dx[1] * u_elem[1] + dN_dx[2] * u_elem[3] +
dN_dx[3] * u_elem[5] + dN_dx[4] * u_elem[7]
εyy = dN_dy[1] * u_elem[2] + dN_dy[2] * u_elem[4] +
dN_dy[3] * u_elem[6] + dN_dy[4] * u_elem[8]
# Engineering shear strain γxy (factor of 2 handled by SymmetricTensor constructor)
γxy = (dN_dy[1] * u_elem[1] + dN_dx[1] * u_elem[2] +
dN_dy[2] * u_elem[3] + dN_dx[2] * u_elem[4] +
dN_dy[3] * u_elem[5] + dN_dx[3] * u_elem[6] +
dN_dy[4] * u_elem[7] + dN_dx[4] * u_elem[8])
# SymmetricTensor{2,2} constructor: (ε11, ε12, ε22)
# Note: ε12 = γxy/2 (tensorial shear strain, not engineering)
return SymmetricTensor{2,2}((εxx, γxy / 2, εyy))
end
@inline function compute_B_transpose_sigma(dN_dx, dN_dy, σ::SymmetricTensor{2,2})
"""
Compute Bᵀ·σ for element residual.
Returns SVector{8} of nodal forces.
"""
# Extract stress components
σxx = σ[1, 1]
σyy = σ[2, 2]
σxy = σ[1, 2] # Tensorial (symmetric), not engineering
# Bᵀ·σ gives forces at each DOF
r_elem = SA[
dN_dx[1]*σxx+dN_dy[1]*σxy, # Node 1, x-direction
dN_dy[1]*σyy+dN_dx[1]*σxy, # Node 1, y-direction
dN_dx[2]*σxx+dN_dy[2]*σxy, # Node 2, x-direction
dN_dy[2]*σyy+dN_dx[2]*σxy, # Node 2, y-direction
dN_dx[3]*σxx+dN_dy[3]*σxy, # Node 3, x-direction
dN_dy[3]*σyy+dN_dx[3]*σxy, # Node 3, y-direction
dN_dx[4]*σxx+dN_dy[4]*σxy, # Node 4, x-direction
dN_dy[4]*σyy+dN_dx[4]*σxy # Node 4, y-direction
]
return r_elem
end
# Main GPU kernel
function elasticity_residual_kernel_tensors!(
r_global::CuDeviceVector{T},
u_global::CuDeviceVector{T},
elem_nodes::CuDeviceMatrix{Int32},
node_coords::CuDeviceMatrix{T},
E::T,
ν::T
) where T
"""
Compute residual using Tensors.jl for proper tensor operations.
Key: ε and σ are SymmetricTensor{2,2}, not vectors!
"""
elem_id = threadIdx().x + (blockIdx().x - 1) * blockDim().x
if elem_id > size(elem_nodes, 1)
return
end
# Material model
material = LinearElastic(E, ν)
# Get element nodes
n1, n2, n3, n4 = elem_nodes[elem_id, 1], elem_nodes[elem_id, 2],
elem_nodes[elem_id, 3], elem_nodes[elem_id, 4]
# Node coordinates
x_coords = SA[node_coords[n1, 1], node_coords[n2, 1],
node_coords[n3, 1], node_coords[n4, 1]]
y_coords = SA[node_coords[n1, 2], node_coords[n2, 2],
node_coords[n3, 2], node_coords[n4, 2]]
# Element DOFs
u_elem = SA[
u_global[2*n1-1], u_global[2*n1],
u_global[2*n2-1], u_global[2*n2],
u_global[2*n3-1], u_global[2*n3],
u_global[2*n4-1], u_global[2*n4]
]
# Accumulate element residual
r_elem = MVector{8,T}(zeros(8))
# Integration loop
for ip in 1:4
ξ, η = GAUSS_POINTS_2D[ip]
w = GAUSS_WEIGHTS_2D[ip]
# Shape function derivatives
dN_dξ, dN_dη = shape_derivatives_quad4(ξ, η)
# Jacobian
J = compute_jacobian_quad4(dN_dξ, dN_dη, x_coords, y_coords)
det_J = det(J)
inv_J = inv(J)
# Physical derivatives: dN/dx = inv(J) · dN/dξ
dN_dx = SA[
inv_J[1, 1]*dN_dξ[1]+inv_J[1, 2]*dN_dη[1],
inv_J[1, 1]*dN_dξ[2]+inv_J[1, 2]*dN_dη[2],
inv_J[1, 1]*dN_dξ[3]+inv_J[1, 2]*dN_dη[3],
inv_J[1, 1]*dN_dξ[4]+inv_J[1, 2]*dN_dη[4]
]
dN_dy = SA[
inv_J[2, 1]*dN_dξ[1]+inv_J[2, 2]*dN_dη[1],
inv_J[2, 1]*dN_dξ[2]+inv_J[2, 2]*dN_dη[2],
inv_J[2, 1]*dN_dξ[3]+inv_J[2, 2]*dN_dη[3],
inv_J[2, 1]*dN_dξ[4]+inv_J[2, 2]*dN_dη[4]
]
# ✅ Compute strain as SymmetricTensor{2,2}
ε = compute_B_matrix_strain(dN_dx, dN_dy, u_elem)
# ✅ Compute stress using material model
σ = compute_stress_2d(material, ε)
# ✅ Compute Bᵀ·σ (element forces)
r_contrib = compute_B_transpose_sigma(dN_dx, dN_dy, σ)
# Accumulate
r_elem .+= r_contrib .* (w * det_J)
end
# Atomic scatter
CUDA.@atomic r_global[2*n1-1] += r_elem[1]
CUDA.@atomic r_global[2*n1] += r_elem[2]
CUDA.@atomic r_global[2*n2-1] += r_elem[3]
CUDA.@atomic r_global[2*n2] += r_elem[4]
CUDA.@atomic r_global[2*n3-1] += r_elem[5]
CUDA.@atomic r_global[2*n3] += r_elem[6]
CUDA.@atomic r_global[2*n4-1] += r_elem[7]
CUDA.@atomic r_global[2*n4] += r_elem[8]
return nothing
end
# ============================================================================
# GPU Assembly Functions (same as before)
# ============================================================================
function compute_residual_gpu!(
r_gpu::CuVector{T},
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T
) where T
n_elements = size(elem_nodes_gpu, 1)
threads = 256
blocks = cld(n_elements, threads)
fill!(r_gpu, zero(T))
@cuda threads = threads blocks = blocks elasticity_residual_kernel_tensors!(
r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν
)
CUDA.synchronize()
return nothing
end
function compute_Jv_gpu!(
Jv_gpu::CuVector{T},
u_gpu::CuVector{T},
v_gpu::CuVector{T},
r0_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
ε::T=T(1e-7)
) where T
u_perturbed = u_gpu .+ ε .* v_gpu
r_perturbed = CUDA.zeros(T, length(u_gpu))
compute_residual_gpu!(r_perturbed, u_perturbed, elem_nodes_gpu, coords_gpu, E, ν)
Jv_gpu .= (r_perturbed .- r0_gpu) ./ ε
return nothing
end
# ============================================================================
# Matrix-Free Operator
# ============================================================================
struct GPUMatrixFreeOperator{T}
u::CuVector{T}
r0::CuVector{T}
elem_nodes::CuMatrix{Int32}
coords::CuMatrix{T}
E::T
ν::T
n::Int
end
Base.size(op::GPUMatrixFreeOperator) = (op.n, op.n)
function LinearAlgebra.mul!(Jv, op::GPUMatrixFreeOperator{T}, v) where T
v_gpu = CuVector{T}(v)
Jv_gpu = CuVector{T}(undef, length(v))
compute_Jv_gpu!(Jv_gpu, op.u, v_gpu, op.r0, op.elem_nodes, op.coords, op.E, op.ν)
copyto!(Jv, Array(Jv_gpu))
return Jv
end
# ============================================================================
# Newton-Krylov Solver
# ============================================================================
function solve_newton_krylov_gpu!(
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
fixed_dofs::Vector{Int};
max_iter::Int=20,
tol::T=T(1e-8),
gmres_tol::T=T(1e-6),
verbose::Bool=true
) where T
n_dofs = length(u_gpu)
r_gpu = CUDA.zeros(T, n_dofs)
for iter in 1:max_iter
compute_residual_gpu!(r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν)
# Enforce BC
r_cpu_temp = Array(r_gpu)
r_cpu_temp[fixed_dofs] .= 0.0
copyto!(r_gpu, r_cpu_temp)
r_norm = CUDA.norm(r_gpu)
if verbose
println(" Newton iter $iter: ||r|| = $r_norm")
end
if r_norm < tol
if verbose
println(" ✅ Converged in $iter iterations")
end
return iter
end
op = GPUMatrixFreeOperator(u_gpu, r_gpu, elem_nodes_gpu, coords_gpu, E, ν, n_dofs)
r_cpu = Array(-r_gpu)
du_cpu, stats = gmres(op, r_cpu, atol=gmres_tol, rtol=0.0, verbose=0)
if !stats.solved
@warn "GMRES did not converge at iteration $iter"
end
du_cpu[fixed_dofs] .= 0.0
du_gpu = CuVector{T}(du_cpu)
u_gpu .+= du_gpu
end
@warn "Newton did not converge in $max_iter iterations"
return max_iter
end
# ============================================================================
# Main Demo
# ============================================================================
function main()
println("\n" * "="^70)
println("GPU Assembly POC with Tensors.jl")
println("="^70)
nx, ny = 10, 10
Lx, Ly = 1.0, 1.0
E, ν = 200e9, 0.3
println("\n📐 Mesh: $(nx*ny) Quad4 elements, $((nx+1)*(ny+1)) nodes, $(2*(nx+1)*(ny+1)) DOFs")
println("🔧 Material: E=$(E/1e9) GPa, ν=$ν (LinearElastic)")
println("✅ Using Tensors.jl: SymmetricTensor{2,2} for ε and σ")
coords, connectivity = generate_rectangle_mesh(nx, ny, Lx, Ly)
n_dofs = 2 * size(coords, 1)
# Boundary conditions
fixed_dofs = Int[]
for node_id in 1:size(coords, 1)
if coords[node_id, 1] < 1e-10
push!(fixed_dofs, 2 * node_id - 1, 2 * node_id)
end
end
u0 = randn(n_dofs) * 1e-6
u0[fixed_dofs] .= 0.0
for node_id in 1:size(coords, 1)
if abs(coords[node_id, 1] - Lx) < 1e-10
u0[2*node_id-1] = 0.001 # 1mm tension
end
end
println("🔒 BC: $(length(fixed_dofs)) fixed DOFs, 1mm tension on right edge")
# Transfer to GPU
elem_nodes_gpu = CuArray{Int32}(connectivity)
coords_gpu = CuArray{Float64}(coords)
u_gpu = CuArray{Float64}(u0)
println("\n🚀 Starting GPU Newton-Krylov (Tensors.jl version)...")
n_iter = solve_newton_krylov_gpu!(
u_gpu, elem_nodes_gpu, coords_gpu, E, ν, fixed_dofs,
max_iter=20, tol=1e-8, gmres_tol=1e-6, verbose=true
)
u_final = Array(u_gpu)
println("\n📊 Results:")
println(" Iterations: $n_iter")
println(" ||u||: $(norm(u_final))")
println("\n✅ POC COMPLETE - Now using proper Tensors.jl!")
println("="^70 * "\n")
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
-430
View File
@@ -1,430 +0,0 @@
"""
GPU Assembly for Tet10 (Quadratic Tetrahedron)
==============================================
The workhorse element for 3D real simulations!
Following JuliaFEM's established pattern:
- Loop through shape functions (i = 1:nnodes)
- Fill 3x3 blocks directly from derivatives dN[1:3, i]
- No "B-matrix" concept - just derivatives!
- Tensors.jl for strain/stress (SymmetricTensor{2,3})
"""
using CUDA
using LinearAlgebra
using Tensors
using Krylov
# ============================================================================
# Material Model
# ============================================================================
struct LinearElastic
E::Float64
ν::Float64
end
@inline λ(mat::LinearElastic) = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
@inline μ(mat::LinearElastic) = mat.E / (2(1 + mat.ν))
"""
3D linear elastic stress: sigma = lambda*tr(eps)*I + 2*mu*eps
"""
@inline function compute_stress_3d(
material::LinearElastic,
eps::SymmetricTensor{2,3,T}
) where T
lambda_val = T(λ(material))
mu_val = T(μ(material))
I = one(eps)
sigma = lambda_val * tr(eps) * I + 2 * mu_val * eps
return sigma
end
# ============================================================================
# Tet10 Reference Element (Quadratic Tetrahedron)
# ============================================================================
"""
Tet10 node numbering (ABAQUS convention):
Vertices: 1-4
Edge midpoints: 5(1-2), 6(2-3), 7(1-3), 8(1-4), 9(2-4), 10(3-4)
"""
const TET10_REF_COORDS = (
Vec{3}((0.0, 0.0, 0.0)), # 1
Vec{3}((1.0, 0.0, 0.0)), # 2
Vec{3}((0.0, 1.0, 0.0)), # 3
Vec{3}((0.0, 0.0, 1.0)), # 4
Vec{3}((0.5, 0.0, 0.0)), # 5 (1-2)
Vec{3}((0.5, 0.5, 0.0)), # 6 (2-3)
Vec{3}((0.0, 0.5, 0.0)), # 7 (1-3)
Vec{3}((0.0, 0.0, 0.5)), # 8 (1-4)
Vec{3}((0.5, 0.0, 0.5)), # 9 (2-4)
Vec{3}((0.0, 0.5, 0.5)) # 10 (3-4)
)
"""
Gauss quadrature for Tet10: 4-point scheme
"""
const GAUSS_TET4 = (
(Vec{3}((0.5854101966249685, 0.1381966011250105, 0.1381966011250105)), 0.25),
(Vec{3}((0.1381966011250105, 0.5854101966249685, 0.1381966011250105)), 0.25),
(Vec{3}((0.1381966011250105, 0.1381966011250105, 0.5854101966249685)), 0.25),
(Vec{3}((0.1381966011250105, 0.1381966011250105, 0.1381966011250105)), 0.25)
)
@inline function tet10_shape_derivatives(xi, eta, zeta)
"""
Quadratic shape function derivatives for Tet10.
Returns tuple of 10 Vec{3} (using Tensors.jl).
"""
lambda = 1 - xi - eta - zeta
# Vertex nodes (1-4)
dN1 = Vec{3}((4 * lambda - 1, 4 * lambda - 1, 4 * lambda - 1))
dN2 = Vec{3}((4 * xi - 1, 0.0, 0.0))
dN3 = Vec{3}((0.0, 4 * eta - 1, 0.0))
dN4 = Vec{3}((0.0, 0.0, 4 * zeta - 1))
# Edge midpoints (5-10)
dN5 = Vec{3}((4 * (1 - 2 * xi - eta - zeta), -4 * xi, -4 * xi))
dN6 = Vec{3}((4 * eta, 4 * xi, 0.0))
dN7 = Vec{3}((-4 * eta, 4 * (1 - xi - 2 * eta - zeta), -4 * eta))
dN8 = Vec{3}((-4 * zeta, -4 * zeta, 4 * (1 - xi - eta - 2 * zeta)))
dN9 = Vec{3}((4 * zeta, 0.0, 4 * xi))
dN10 = Vec{3}((0.0, 4 * zeta, 4 * eta))
return (dN1, dN2, dN3, dN4, dN5, dN6, dN7, dN8, dN9, dN10)
end
@inline function compute_jacobian_tet10(dN_dxi::NTuple{10,Vec{3,T}}, X::NTuple{10,Vec{3,T}}) where T
"""
Jacobian using tensor products: J = Σ_i dN_i ⊗ X_i
"""
return sum(dN_dxi[i] X[i] for i in 1:10)
end
# ============================================================================
# JuliaFEM Pattern: Loop Through Shape Functions, Fill 3x3 Blocks
# ============================================================================
@inline function compute_strain_from_displacements(
dN_dx::NTuple{10,Vec{3,T}},
u::NTuple{10,Vec{3,T}}
) where T
"""
Compute strain using Tensors.jl tensor products.
Displacement gradient: ∇u = Σ_i dN_i ⊗ u_i
Strain (small): ε = 1/2 (∇u + ∇uᵀ) = sym(∇u)
"""
# Displacement gradient via tensor products
gradu = sum(dN_dx[i] u[i] for i in 1:10)
# Symmetric part (strain)
return symmetric(gradu)
end
@inline function compute_nodal_forces_from_stress(
dN_dx::NTuple{10,Vec{3,T}},
sigma::SymmetricTensor{2,3,T}
) where T
"""
Compute nodal forces using Tensors.jl.
Force at node i: f_i = dN_i · σ
"""
return ntuple(i -> dN_dx[i] sigma, Val(10))
end
# ============================================================================
# GPU Kernel
# ============================================================================
function tet10_residual_kernel!(
r_global::CuDeviceVector{T},
u_global::CuDeviceVector{T},
elem_nodes::CuDeviceMatrix{Int32},
node_coords::CuDeviceMatrix{T},
E::T,
ν::T
) where T
"""
Element-parallel GPU kernel for Tet10 elasticity.
Pattern:
1. Get 10 node coordinates and 30 DOFs
2. Loop over 4 Gauss points
3. Compute shape derivatives dN/dξ
4. Compute Jacobian and physical derivatives dN/dx
5. Compute strain from derivatives (no B-matrix!)
6. Compute stress from material model
7. Compute nodal forces from stress (loop through shape functions)
8. Atomic scatter to global residual
"""
elem_id = threadIdx().x + (blockIdx().x - 1) * blockDim().x
if elem_id > size(elem_nodes, 1)
return
end
material = LinearElastic(E, ν)
# Get element nodes (1-indexed to 10 nodes)
nodes = ntuple(i -> elem_nodes[elem_id, i], Val(10))
# Element coordinates (10 nodes × 3 coordinates)
X = ntuple(Val(10)) do i
Vec{3}((node_coords[nodes[i], 1],
node_coords[nodes[i], 2],
node_coords[nodes[i], 3]))
end
# Element displacements (10 nodes × 3 DOFs)
u = ntuple(Val(10)) do i
Vec{3}((u_global[3*nodes[i]-2],
u_global[3*nodes[i]-1],
u_global[3*nodes[i]]))
end
# Accumulate element residual (10 forces as Vec{3})
r_elem = [zero(Vec{3,T}) for _ in 1:10]
# Integration loop (4 Gauss points)
for (xez_tuple, w) in GAUSS_TET4
xi, eta, zeta = xez_tuple[1], xez_tuple[2], xez_tuple[3]
# Shape function derivatives in reference coordinates
dN_dxi = tet10_shape_derivatives(xi, eta, zeta)
# Jacobian using tensor products: J = Σ dN_i ⊗ X_i
J = compute_jacobian_tet10(dN_dxi, X)
detJ = det(J)
invJ = inv(J)
# Physical derivatives: dN_dx = invJ · dN_dxi (tensor contraction)
dN_dx = ntuple(i -> invJ dN_dxi[i], Val(10))
# Compute strain: ε = sym(∇u) where ∇u = Σ dN_i ⊗ u_i
eps = compute_strain_from_displacements(dN_dx, u)
# Compute stress
sigma = compute_stress_3d(material, eps)
# Compute nodal forces: f_i = dN_i · σ
f_contrib = compute_nodal_forces_from_stress(dN_dx, sigma)
# Accumulate with quadrature weight
for i in 1:10
r_elem[i] += f_contrib[i] * (w * detJ)
end
end
# Atomic scatter (10 nodes × 3 components)
for i in 1:10
CUDA.@atomic r_global[3*nodes[i]-2] += r_elem[i][1]
CUDA.@atomic r_global[3*nodes[i]-1] += r_elem[i][2]
CUDA.@atomic r_global[3*nodes[i]] += r_elem[i][3]
end
return nothing
end
# ============================================================================
# Assembly & Solver Wrappers
# ============================================================================
function compute_residual_gpu!(
r_gpu::CuVector{T},
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T
) where T
n_elements = size(elem_nodes_gpu, 1)
threads = 256
blocks = cld(n_elements, threads)
fill!(r_gpu, zero(T))
@cuda threads = threads blocks = blocks tet10_residual_kernel!(
r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν
)
CUDA.synchronize()
return nothing
end
function compute_Jv_gpu!(
Jv_gpu::CuVector{T},
u_gpu::CuVector{T},
v_gpu::CuVector{T},
r0_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
ε::T=T(1e-7)
) where T
u_perturbed = u_gpu .+ ε .* v_gpu
r_perturbed = CUDA.zeros(T, length(u_gpu))
compute_residual_gpu!(r_perturbed, u_perturbed, elem_nodes_gpu, coords_gpu, E, ν)
Jv_gpu .= (r_perturbed .- r0_gpu) ./ ε
return nothing
end
struct GPUMatrixFreeOperator{T}
u::CuVector{T}
r0::CuVector{T}
elem_nodes::CuMatrix{Int32}
coords::CuMatrix{T}
E::T
ν::T
n::Int
end
Base.size(op::GPUMatrixFreeOperator) = (op.n, op.n)
function LinearAlgebra.mul!(Jv, op::GPUMatrixFreeOperator{T}, v) where T
v_gpu = CuVector{T}(v)
Jv_gpu = CuVector{T}(undef, length(v))
compute_Jv_gpu!(Jv_gpu, op.u, v_gpu, op.r0, op.elem_nodes, op.coords, op.E, op.ν)
copyto!(Jv, Array(Jv_gpu))
return Jv
end
function solve_newton_krylov_gpu!(
u_gpu::CuVector{T},
elem_nodes_gpu::CuMatrix{Int32},
coords_gpu::CuMatrix{T},
E::T,
ν::T,
fixed_dofs::Vector{Int};
max_iter::Int=20,
tol::T=T(1e-8),
gmres_tol::T=T(1e-6),
verbose::Bool=true
) where T
n_dofs = length(u_gpu)
r_gpu = CUDA.zeros(T, n_dofs)
for iter in 1:max_iter
compute_residual_gpu!(r_gpu, u_gpu, elem_nodes_gpu, coords_gpu, E, ν)
# Enforce BC
r_cpu_temp = Array(r_gpu)
r_cpu_temp[fixed_dofs] .= 0.0
copyto!(r_gpu, r_cpu_temp)
r_norm = CUDA.norm(r_gpu)
if verbose
println(" Newton iter $iter: ||r|| = $r_norm")
end
if r_norm < tol
if verbose
println(" ✅ Converged in $iter iterations")
end
return iter
end
op = GPUMatrixFreeOperator(u_gpu, r_gpu, elem_nodes_gpu, coords_gpu, E, ν, n_dofs)
r_cpu = Array(-r_gpu)
du_cpu, stats = gmres(op, r_cpu, atol=gmres_tol, rtol=0.0, verbose=0)
if !stats.solved
@warn "GMRES did not converge at iteration $iter"
end
du_cpu[fixed_dofs] .= 0.0
du_gpu = CuVector{T}(du_cpu)
u_gpu .+= du_gpu
end
@warn "Newton did not converge in $max_iter iterations"
return max_iter
end
# ============================================================================
# Test Mesh Generation
# ============================================================================
function generate_single_tet10_mesh()
"""
Single Tet10 element for testing.
"""
# 10 nodes: 4 vertices + 6 edge midpoints
coords = [
0.0 0.0 0.0;
1.0 0.0 0.0;
0.0 1.0 0.0;
0.0 0.0 1.0;
0.5 0.0 0.0;
0.5 0.5 0.0;
0.0 0.5 0.0;
0.0 0.0 0.5;
0.5 0.0 0.5;
0.0 0.5 0.5
]
connectivity = reshape(Int32[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 1, 10)
return coords, connectivity
end
# ============================================================================
# Main Demo
# ============================================================================
function main()
println("\n" * "="^70)
println("GPU Assembly: Tet10 (Quadratic Tetrahedron)")
println("="^70)
E, ν = 200e9, 0.3
coords, connectivity = generate_single_tet10_mesh()
n_dofs = 3 * size(coords, 1)
println("\n📐 Mesh: 1 Tet10 element, 10 nodes, 30 DOFs")
println("🔧 Material: E=$(E/1e9) GPa, ν=$ν")
println("✅ Pattern: Loop through shape functions, fill 3x3 blocks")
println("✅ No B-matrix concept - just derivatives!")
# Boundary conditions: Fix node 1 (DOFs 1,2,3)
fixed_dofs = [1, 2, 3]
# Initial displacement: Small perturbation
u0 = randn(n_dofs) * 1e-6
u0[fixed_dofs] .= 0.0
u0[4] = 0.001 # Pull node 2 in x-direction
println("🔒 BC: Node 1 fixed, node 2 displaced 1mm in x")
# Transfer to GPU
elem_nodes_gpu = CuArray{Int32}(connectivity)
coords_gpu = CuArray{Float64}(coords)
u_gpu = CuArray{Float64}(u0)
println("\n🚀 Starting GPU Newton-Krylov (Tet10)...")
n_iter = solve_newton_krylov_gpu!(
u_gpu, elem_nodes_gpu, coords_gpu, E, ν, fixed_dofs,
max_iter=20, tol=1e-8, gmres_tol=1e-6, verbose=true
)
u_final = Array(u_gpu)
println("\n📊 Results:")
println(" Iterations: $n_iter")
println(" ||u||: $(norm(u_final))")
println("\n✅ TET10 POC COMPLETE!")
println("="^70 * "\n")
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
-518
View File
@@ -1,518 +0,0 @@
#!/usr/bin/env julia
#
# GPU ElementSet Demo: Immutable Fields + Type Stability
#
# Demonstrates:
# 1. ElementSet pattern with type-stable fields
# 2. GPU kernel execution with immutable field access
# 3. GENERAL approach: Elements go to GPU with connectivity inside
# 4. Zero-allocation assembly loop pattern
# 5. Creating new field containers between time steps (cheap!)
#
# Key architectural decision:
# - Element struct contains connectivity (NTuple - immutable, zero-cost)
# - Fields live in ElementSet (type-stable, immutable)
# - GPU kernel receives Vector{Element} directly (general!)
# - Connectivity accessed via element.connectivity on GPU
#
# This uses MOCK GPU execution (no CUDA.jl dependency) to show the pattern.
# Real GPU code would look identical - that's the point!
using BenchmarkTools
using LinearAlgebra
using Printf
println("="^70)
println("GPU ElementSet Demo: Immutable Fields Pattern")
println("="^70)
println()
# ============================================================================
# Mock Structures (Simplified for Demonstration)
# ============================================================================
"""Mock element - simple, no fields inside"""
struct Element{N,B}
id::UInt
connectivity::NTuple{N,UInt}
basis::B
end
"""Mock basis type"""
struct MockBasis end
"""Element set - elements + type-stable fields"""
struct ElementSet{E,F}
name::String
elements::Vector{E}
fields::F # Type-stable! Can be NamedTuple, struct, anything
end
"""Assembly cache - pre-allocated buffers for zero-allocation assembly"""
struct AssemblyCache
K_local::Matrix{Float64}
f_local::Vector{Float64}
N::Vector{Float64}
dN::Matrix{Float64}
end
function AssemblyCache(ndof_local::Int)
return AssemblyCache(
zeros(ndof_local, ndof_local),
zeros(ndof_local),
zeros(4), # 4 basis functions for quad
zeros(2, 4), # 2D derivatives
)
end
# ============================================================================
# Mock GPU Module (Simulates CUDA.jl behavior)
# ============================================================================
module MockGPU
"""Mock GPU array type"""
struct CuArray{T,N}
data::Array{T,N}
end
Base.length(a::CuArray) = length(a.data)
Base.getindex(a::CuArray, i...) = getindex(a.data, i...)
Base.setindex!(a::CuArray, v, i...) = setindex!(a.data, v, i...)
"""Transfer to GPU (mock - just wraps array)"""
cu(x::Array) = CuArray(x)
cu(x::Vector) = CuArray(x) # Also handle vectors
"""Transfer from GPU (mock - unwraps array)"""
cpu(x::CuArray) = x.data
"""Mock @cuda macro - just calls function serially"""
macro cuda(args...)
# Extract function call from threads=... blocks=... function(...)
func_call = args[end]
return esc(quote
# In real CUDA, this would launch kernel on GPU
# Here we just call it serially to show the pattern
$func_call
end)
end
export CuArray, cu, cpu, @cuda
end
using .MockGPU# ============================================================================
# GPU Kernel: Element Assembly (Type-Stable!)
# ============================================================================
"""
GPU kernel for element assembly - GENERAL VERSION.
Works with Element struct directly (connectivity inside elements).
Key points:
1. Elements vector is transferred to GPU (struct-of-arrays pattern)
2. Fields are immutable (read-only access)
3. Type-stable: Element{N,B} has known connectivity length N
4. Zero allocations (connectivity is NTuple, immutable)
In real CUDA: Each thread processes one element
"""
function gpu_assemble_kernel!(
K_global::CuArray{Float64,2},
f_global::CuArray{Float64,1},
elements::CuArray{Element{4,MockBasis},1}, # Vector of elements
E::Float64, # Young's modulus (immutable constant)
ν::Float64, # Poisson's ratio (immutable constant)
u::CuArray{Float64,2}, # Displacement (immutable, read-only)
n_elements::Int,
)
# In real CUDA: thread_id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
# Here we loop serially to simulate
for elem_id in 1:n_elements
# Get element from GPU array (elements were transferred!)
element = elements[elem_id]
# Get connectivity from element struct (NTuple - zero allocation!)
# This is the key: connectivity lives INSIDE the element
nodes = element.connectivity
# Mock assembly computation
# In real code: compute K_local from E, ν, u[nodes]
# Here we just do simple arithmetic to show the pattern
K_local_value = E * (1 - ν^2) # Mock stiffness
# Mock: add contribution to diagonal (in real code: full K_local)
for i in 1:4
node = nodes[i]
# Atomic add in real CUDA
K_global[node, node] += K_local_value * 0.25
f_global[node] += K_local_value * u[1, node] * 0.1
end
end
return nothing
end
# ============================================================================
# CPU Assembly (Same Logic, No GPU)
# ============================================================================
"""CPU version of assembly - same logic as GPU kernel"""
function cpu_assemble!(
K_global::Matrix{Float64},
f_global::Vector{Float64},
element_set::ElementSet,
cache::AssemblyCache,
)
# Access fields (type-stable!)
fields = element_set.fields
E = fields.E
ν = fields.ν
u = fields.u
# Assembly loop (zero allocations!)
for element in element_set.elements
nodes = element.connectivity
# Mock assembly
K_local_value = E * (1 - ν^2)
for i in 1:length(nodes)
node = nodes[i]
K_global[node, node] += K_local_value * 0.25
f_global[node] += K_local_value * u[1, node] * 0.1
end
end
return nothing
end
# ============================================================================
# Setup Problem
# ============================================================================
println("Setting up problem...")
println()
# Problem size
n_nodes = 1000
n_elements = 800 # Quad elements, 4 nodes each
# Create elements (no fields inside!)
elements = [
Element{4,MockBasis}(
UInt(i),
(UInt(i), UInt(i + 1), UInt(i + 101), UInt(i + 100)), # Mock connectivity
MockBasis()
)
for i in 1:n_elements
]
# Initial fields (immutable!)
fields_initial = (
E=210e3, # Young's modulus
ν=0.3, # Poisson's ratio
u=zeros(3, n_nodes), # Initial displacement
)
# Create element set
element_set = ElementSet("steel_body", elements, fields_initial)
println("Problem setup:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" Element type: ", typeof(elements[1]))
println(" Field type: ", typeof(element_set.fields))
println(" Fields are immutable: ", !ismutable(element_set.fields))
println()
# ============================================================================
# CPU Assembly Benchmark
# ============================================================================
println("="^70)
println("CPU Assembly (Baseline)")
println("="^70)
println()
K_cpu = zeros(n_nodes, n_nodes)
f_cpu = zeros(n_nodes)
cache = AssemblyCache(12) # 3 DOF × 4 nodes
println("First assembly (with compilation):")
@time cpu_assemble!(K_cpu, f_cpu, element_set, cache)
println("\nBenchmarked assembly:")
cpu_result = @benchmark cpu_assemble!($K_cpu, $f_cpu, $element_set, $cache) setup = (
K_cpu = zeros($n_nodes, $n_nodes);
f_cpu = zeros($n_nodes)
)
display(cpu_result)
println()
cpu_time = median(cpu_result).time / 1e6 # Convert to ms
cpu_allocs = median(cpu_result).allocs
println("\n📊 CPU Results:")
@printf(" Time: %.3f ms\n", cpu_time)
println(" Allocations: $cpu_allocs")
println(" ✓ Zero allocations in assembly loop: ", cpu_allocs == 0 ? "YES ✅" : "NO ❌")
println()
# ============================================================================
# GPU Assembly (Mock)
# ============================================================================
println("="^70)
println("GPU Assembly (Mock CUDA)")
println("="^70)
println()
# Transfer data to GPU
println("Transferring data to GPU...")
println(" Key insight: Elements themselves go to GPU!")
println(" Connectivity lives INSIDE each element (NTuple)")
println()
# Transfer elements to GPU (GENERAL APPROACH!)
elements_gpu = cu(elements) # Vector{Element{4,MockBasis}} → CuArray
# Transfer field data
u_gpu = cu(fields_initial.u)
K_gpu = cu(zeros(n_nodes, n_nodes))
f_gpu = cu(zeros(n_nodes))
println(" Elements: ", typeof(elements_gpu))
println(" Displacement: ", typeof(u_gpu))
println(" Stiffness: ", typeof(K_gpu))
println()
# Launch kernel (GENERAL VERSION - takes elements directly!)
println("Launching GPU kernel...")
@cuda threads = 256 blocks = ceil(Int, n_elements / 256) gpu_assemble_kernel!(
K_gpu, f_gpu, elements_gpu, # ← Elements, not separate connectivity!
fields_initial.E, fields_initial.ν, u_gpu,
n_elements
)
println(" ✓ Kernel execution complete")
println(" ✓ Elements accessed directly on GPU")
println(" ✓ Connectivity read from element.connectivity")
println()
# Transfer back
K_gpu_result = cpu(K_gpu)
f_gpu_result = cpu(f_gpu)
# Verify correctness
K_cpu_check = zeros(n_nodes, n_nodes)
f_cpu_check = zeros(n_nodes)
cpu_assemble!(K_cpu_check, f_cpu_check, element_set, cache)
error_K = norm(K_gpu_result - K_cpu_check) / (norm(K_cpu_check) + 1e-10)
error_f = norm(f_gpu_result - f_cpu_check) / (norm(f_cpu_check) + 1e-10)
println("📊 GPU Results:")
println(" Relative error (K): ", @sprintf("%.6e", error_K))
println(" Relative error (f): ", @sprintf("%.6e", error_f))
println(" ✓ GPU matches CPU: ", (error_K < 1e-6 && error_f < 1e-6) ? "YES ✅" : "NO ❌")
println()
# ============================================================================
# Time Stepping: Creating New Field Containers
# ============================================================================
println("="^70)
println("Time Stepping: Immutable Fields Pattern")
println("="^70)
println()
println("Simulating 5 time steps with field updates...")
println()
# Simulate some displacement changes
displacement_changes = [
0.001 * sin(2π * t) * ones(3, n_nodes) for t in range(0, 1, length=5)
]
time_step_times = Float64[]
# Need to track element_set explicitly for time stepping
global current_element_set = element_set
for (step, du) in enumerate(displacement_changes)
println("Time step $step:")
# Get current displacement
u_old = current_element_set.fields.u
# Compute new displacement (mock solver)
u_new = u_old + du
# Create NEW field container (immutable pattern!)
# This is CHEAP - just wraps references, no copying!
fields_new = (
E=current_element_set.fields.E, # Keep old (constant)
ν=current_element_set.fields.ν, # Keep old (constant)
u=u_new, # New displacement
)
# Create new element set (also cheap - just wraps references)
element_set_new = ElementSet(current_element_set.name, current_element_set.elements, fields_new)
# Assemble with new fields
K_new = zeros(n_nodes, n_nodes)
f_new = zeros(n_nodes)
time_step = @elapsed cpu_assemble!(K_new, f_new, element_set_new, cache)
push!(time_step_times, time_step * 1000) # Convert to ms
@printf(" Assembly time: %.3f ms\n", time_step * 1000)
println(" Max displacement: ", @sprintf("%.6e", maximum(abs.(u_new))))
println(" Field container recreated: ✓")
println()
# Update for next iteration
global current_element_set = element_set_new
end
println("📊 Time Stepping Results:")
@printf(" Average assembly time: %.3f ms\n", sum(time_step_times) / length(time_step_times))
println(" Field updates: ", length(displacement_changes))
println(" ✓ Creating new field containers is cheap (no copying)")
println()
# ============================================================================
# Memory Usage Analysis
# ============================================================================
println("="^70)
println("Memory Usage: Immutable vs Mutable")
println("="^70)
println()
# Immutable pattern
fields_immutable = (E=210e3, ν=0.3, u=zeros(3, n_nodes))
size_immutable = sizeof(fields_immutable) + sizeof(fields_immutable.u)
# Hypothetical mutable pattern
mutable struct MutableFields
E::Float64
ν::Float64
u::Matrix{Float64}
end
fields_mutable = MutableFields(210e3, 0.3, zeros(3, n_nodes))
size_mutable = sizeof(fields_mutable) + sizeof(fields_mutable.u)
println("Memory comparison:")
@printf(" Immutable (NamedTuple): %d bytes\n", size_immutable)
@printf(" Mutable (struct): %d bytes\n", size_mutable)
@printf(" Difference: %.1f%%\n", 100 * (size_immutable - size_mutable) / size_mutable)
println()
println("Creating new containers (benchmark):")
fields_base = (E=210e3, ν=0.3, u=zeros(3, n_nodes))
u_sample = zeros(3, n_nodes)
println("\nImmutable pattern (create new NamedTuple):")
@btime (E=$fields_base.E, ν=$fields_base.ν, u=$u_sample)
println("\nMutable pattern (would need to copy for safety):")
@btime deepcopy($fields_mutable)
println()
println("Key insight: Creating new NamedTuple is ~1000× faster than deepcopy!")
println(" (NamedTuple just wraps references, no data copying)")
println()
# ============================================================================
# Summary and Validation
# ============================================================================
println("="^70)
println("SUMMARY: ElementSet + Immutable Fields Pattern")
println("="^70)
println()
validation_passed = true
println("✓ Element structure:")
println(" - No fields inside Element struct")
println(" - Connectivity is NTuple (zero-cost)")
println(" - Simple, type-stable")
println()
println("✓ ElementSet structure:")
println(" - Groups elements + fields")
println(" - Fields type parameter F is type-stable")
println(" - Works with NamedTuple, custom struct, anything")
println()
println("✓ GPU compatibility:")
if error_K < 1e-10 && error_f < 1e-10
println(" - GPU kernel executed successfully ✅")
println(" - Results match CPU (error < 1e-10)")
println(" - Immutable field access works on GPU")
println(" - GENERAL: Elements transferred directly to GPU")
println(" - Connectivity accessed from element.connectivity")
println(" - No manual data extraction needed!")
else
println(" - GPU execution had errors ❌")
validation_passed = false
end
println()
println("✓ Performance:")
if cpu_allocs == 0
println(" - Zero allocations in assembly loop ✅")
@printf(" - Assembly time: %.3f ms for %d elements\n", cpu_time, n_elements)
else
println(" - Assembly had allocations ❌ ($cpu_allocs)")
validation_passed = false
end
println()
println("✓ Immutability pattern:")
println(" - Fields are immutable (NamedTuple)")
println(" - Create new containers between time steps")
println(" - Creating new NamedTuple: ~10 ns (no copying!)")
println(" - Deepcopy mutable struct: ~10 μs (1000× slower)")
println()
println("✓ Code clarity:")
println(" - Same code for CPU and GPU (type-stable)")
println(" - ElementSet matches physical thinking (properties per set)")
println(" - Separation: Element (geometry) vs Fields (properties)")
println()
println("="^70)
if validation_passed
println("✅ ALL VALIDATIONS PASSED")
println()
println("The ElementSet + immutable fields pattern is:")
println(" 1. Type-stable (9-92× faster than Dict)")
println(" 2. GPU-compatible (works on CUDA)")
println(" 3. Zero-allocation (in assembly loop)")
println(" 4. Physically correct (properties per set)")
println(" 5. Fast to update (creating new containers is cheap)")
else
println("⚠️ SOME VALIDATIONS FAILED")
println("Review implementation details")
end
println("="^70)
println()
println("Next steps:")
println(" 1. Replace this mock with real CUDA.jl")
println(" 2. Implement in src/elements/elements.jl")
println(" 3. Update Problem struct to use ElementSet")
println(" 4. Migrate examples to new pattern")
println(" 5. Add benchmarks to CI")
println()
println("See also:")
println(" - benchmarks/field_storage_comparison.jl (CPU benchmarks)")
println(" - docs/book/element_field_architecture.md (design rationale)")
println(" - demos/gpu_mpi_demo.jl (real multi-GPU example)")
-322
View File
@@ -1,322 +0,0 @@
#!/usr/bin/env julia
#
# GPU and MPI Real Hardware Demonstration
#
# This script demonstrates that type-stable field data flows to:
# 1. Real CUDA GPU - actual GPU kernel execution
# 2. Real MPI processes - actual inter-process communication
#
# Requirements:
# - CUDA-capable GPU (optional, will detect)
# - MPI installation
# - Run with: mpirun -np 2 julia --project=. benchmarks/gpu_mpi_demo.jl
#
# KEY INSIGHT: Type-stable code on CPU translates directly to GPU/MPI.
#
using LinearAlgebra
# Try to load CUDA (optional)
CUDA_AVAILABLE = false
try
using CUDA
if CUDA.functional()
global CUDA_AVAILABLE = true
println("✓ CUDA GPU detected: $(CUDA.name(CUDA.device()))")
else
println("⚠ CUDA.jl installed but no GPU detected")
end
catch e
println(" CUDA.jl not available (optional): $e")
end
# Load MPI (required)
using MPI
MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)
# Only rank 0 prints headers
function println_master(args...)
if rank == 0
println(args...)
end
end
println_master("="^70)
println_master("GPU and MPI Real Hardware Demonstration")
println_master("="^70)
println_master("MPI: rank=$rank/$size")
println_master()
# ============================================================================
# Part 1: Type-Stable Data Structures
# ============================================================================
println_master("Part 1: Type-Stable Data Structures")
println_master("-"^70)
# Define problem data (type-stable!)
n_nodes_per_rank = 1000
n_elements_per_rank = 100
# Each rank owns a subdomain
nodes = rand(Float64, 3, n_nodes_per_rank)
connectivity = rand(1:n_nodes_per_rank, 8, n_elements_per_rank)
E = 210e3 # Young's modulus (Float64)
ν = 0.3 # Poisson's ratio (Float64)
displacement = rand(Float64, 3, n_nodes_per_rank)
println_master("✓ Created typed data structures on each rank:")
println_master(" - nodes: Array{Float64,2}")
println_master(" - connectivity: Array{Int,2}")
println_master(" - displacement: Array{Float64,2}")
println_master(" - E: Float64 = $E")
println_master(" - ν: Float64 = $ν")
println_master()
# ============================================================================
# Part 2: MPI Communication (Real Hardware)
# ============================================================================
println_master("Part 2: MPI Data Transfer Between Ranks")
println_master("-"^70)
# Synchronize all ranks
MPI.Barrier(comm)
if rank == 0
# Rank 0 sends to rank 1
println("Rank 0: Sending displacement data to rank 1...")
nbytes = sizeof(displacement)
MPI.Send(displacement, comm; dest=1, tag=0)
println("Rank 0: Sent $(nbytes) bytes ($(nbytes/1024) KB)")
# Also send material properties
material = [E, ν]
MPI.Send(material, comm; dest=1, tag=1)
println("Rank 0: Sent material properties")
elseif rank == 1 && size >= 2
# Rank 1 receives from rank 0
println("Rank 1: Receiving displacement data from rank 0...")
received_disp = similar(displacement)
MPI.Recv!(received_disp, comm; source=0, tag=0)
nbytes = sizeof(received_disp)
println("Rank 1: Received $(nbytes) bytes ($(nbytes/1024) KB)")
# Receive material properties
received_mat = zeros(Float64, 2)
MPI.Recv!(received_mat, comm; source=0, tag=1)
println("Rank 1: Received material properties: E=$(received_mat[1]), ν=$(received_mat[2])")
# Verify data integrity
checksum = sum(abs, received_disp)
println("Rank 1: Data checksum = $(checksum)")
end
MPI.Barrier(comm)
println_master()
println_master("✅ MPI communication successful!")
println_master(" Type-stable arrays (Matrix{Float64}) transferred efficiently")
println_master()
# ============================================================================
# Part 3: GPU Kernel Execution (Real Hardware, if available)
# ============================================================================
if CUDA_AVAILABLE && rank == 0
println_master("Part 3: GPU Kernel Execution (Real CUDA Hardware)")
println_master("-"^70)
# Define a simple assembly kernel
function assemble_element_kernel!(
K_elements::CuDeviceMatrix{Float64},
nodes::CuDeviceMatrix{Float64},
connectivity::CuDeviceMatrix{Int32},
E::Float64,
ν::Float64,
n_elements::Int32
)
# GPU thread indexing
idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if idx <= n_elements
# Mock stiffness computation
# In real FEM: would access element nodes and integrate
# Here: simplified to show type-stable GPU execution
K_local = E * (1 - ν^2)
# Store result (simplified: single value per element)
K_elements[idx, 1] = K_local
end
return nothing
end
println("✓ Preparing data for GPU transfer:")
println(" - nodes: $(sizeof(nodes)) bytes")
println(" - connectivity: $(sizeof(connectivity)) bytes")
# Transfer to GPU
println("\n✓ Transferring data to GPU...")
d_nodes = CuArray(nodes)
d_connectivity = CuArray(Int32.(connectivity))
d_K_elements = CUDA.zeros(Float64, n_elements_per_rank, 64)
nbytes_transferred = sizeof(nodes) + sizeof(connectivity)
println(" Transferred $(nbytes_transferred) bytes to GPU")
# Launch kernel
println("\n✓ Launching GPU kernel...")
threads_per_block = 256
blocks = cld(n_elements_per_rank, threads_per_block)
@cuda threads = threads_per_block blocks = blocks assemble_element_kernel!(
d_K_elements, d_nodes, d_connectivity, E, ν, Int32(n_elements_per_rank)
)
CUDA.synchronize()
println(" Kernel executed on $(blocks) blocks × $(threads_per_block) threads")
# Transfer results back
println("\n✓ Transferring results from GPU...")
K_elements = Array(d_K_elements)
println(" Transferred $(sizeof(K_elements)) bytes from GPU")
# Verify results
println("\n✓ Verifying results:")
expected_value = E * (1 - ν^2)
actual_value = K_elements[1, 1]
println(" Expected: $(expected_value)")
println(" Actual: $(actual_value)")
println(" Match: $(abs(expected_value - actual_value) < 1e-10 ? "" : "")")
println("\n✅ GPU execution successful!")
println(" Type-stable kernel compiled and executed on real GPU hardware")
println()
elseif rank == 0
println_master("Part 3: GPU Execution")
println_master("-"^70)
println_master(" No CUDA GPU available (optional)")
println_master(" Type-stable code WOULD compile for GPU if hardware present")
println_master()
end
# ============================================================================
# Part 4: Combined GPU + MPI Pattern (if GPU available)
# ============================================================================
if CUDA_AVAILABLE && size >= 2
println_master("Part 4: Combined GPU + MPI Workflow")
println_master("-"^70)
MPI.Barrier(comm)
if rank == 0
println("Rank 0: Computing on GPU...")
# GPU computation
d_result = CUDA.zeros(Float64, n_elements_per_rank)
# (kernel launch would go here)
d_result .= E * (1 - ν^2)
# Transfer back from GPU
cpu_result = Array(d_result)
println("Rank 0: Got results from GPU ($(length(cpu_result)) elements)")
# Send to rank 1 via MPI
println("Rank 0: Sending GPU results to rank 1 via MPI...")
MPI.Send(cpu_result, comm; dest=1, tag=10)
println("Rank 0: Sent $(sizeof(cpu_result)) bytes")
elseif rank == 1
println("Rank 1: Waiting for GPU results from rank 0...")
# Receive from rank 0
received_result = zeros(Float64, n_elements_per_rank)
MPI.Recv!(received_result, comm; source=0, tag=10)
println("Rank 1: Received $(sizeof(received_result)) bytes from rank 0's GPU")
# Verify
checksum = sum(received_result)
println("Rank 1: Result checksum = $(checksum)")
end
MPI.Barrier(comm)
println_master()
println_master("✅ Combined GPU+MPI workflow successful!")
println_master(" Data flowed: Rank 0 GPU → Rank 0 CPU → MPI → Rank 1 CPU")
println_master()
end
# ============================================================================
# Summary
# ============================================================================
MPI.Barrier(comm)
if rank == 0
println("="^70)
println("SUMMARY: Type Stability Enables GPU and MPI")
println("="^70)
println()
println("✅ Demonstrated on Real Hardware:")
println()
println("1. MPI Communication:")
println(" • Transferred Matrix{Float64} between ranks")
println(" • Fast buffer transfer (not serialization)")
println(" • Type: $(typeof(displacement))")
println(" • Size: $(sizeof(displacement)) bytes")
println()
if CUDA_AVAILABLE
println("2. GPU Execution:")
println(" • Compiled type-stable kernel for GPU")
println(" • Executed on real CUDA hardware")
println(" • Zero allocations in kernel")
println(" • Device: $(CUDA.name(CUDA.device()))")
println()
if size >= 2
println("3. Combined Workflow:")
println(" • GPU computation on rank 0")
println(" • MPI transfer to rank 1")
println(" • End-to-end type stability")
println()
end
else
println("2. GPU Execution:")
println(" • No GPU detected (optional)")
println(" • Type-stable code ready for GPU")
println()
end
println("Key Insights:")
println()
println("• Type stability is REQUIRED for GPU compilation")
println(" - Dict{String,Any} would FAIL to compile for GPU")
println(" - Float64, Matrix{Float64} compile successfully")
println()
println("• Type stability enables fast MPI transfers")
println(" - Typed arrays: fast buffer transfer")
println(" - Mixed types: slow serialization (~100× slower)")
println()
println("• Same code pattern works everywhere")
println(" - CPU: 9-92× speedup (measured)")
println(" - GPU: Enables execution (requirement)")
println(" - MPI: Fast transfers (requirement)")
println()
println("CONCLUSION:")
println("Type-stable field storage is not optional—it's the foundation")
println("for high-performance FEM on modern hardware (GPU, MPI, threading).")
println()
println("="^70)
end
MPI.Finalize()
-327
View File
@@ -1,327 +0,0 @@
#!/usr/bin/env julia
#
# GPU and MPI Mock Demonstration
#
# This script demonstrates that type-stable field data can flow to:
# 1. GPU (CUDA) - using mock kernel without requiring CUDA.jl dependency
# 2. MPI processes - showing efficient data transfer patterns
#
# KEY INSIGHT: Type-stable code on CPU translates directly to GPU/MPI.
# The same zero-allocation patterns work across all execution models.
#
println("="^70)
println("GPU and MPI Data Flow Demonstration")
println("="^70)
println()
# ============================================================================
# Mock CUDA Kernel (Minimal CUDA-like interface without dependency)
# ============================================================================
"""
Mock CUDA-like interface demonstrating type-stable kernel compilation.
In real CUDA.jl:
@cuda threads=256 blocks=ceil(Int, n/256) my_kernel!(data, n)
The key requirement: ALL code in kernel must be type-stable.
Type instability (Dict{String,Any}, Any types) causes compilation failure.
"""
module MockCUDA
# Mock CuArray that acts like a typed GPU array
struct CuArray{T,N}
data::Array{T,N} # In reality, this would be device memory
end
# Mock transfer to device
function cu(arr::Array{T,N}) where {T,N}
println(" 📤 Transferring $(sizeof(arr)) bytes to GPU (mock)")
return CuArray{T,N}(copy(arr))
end
# Mock transfer from device
function Array(carr::CuArray{T,N}) where {T,N}
println(" 📥 Transferring $(sizeof(carr.data)) bytes from GPU (mock)")
return copy(carr.data)
end
# Mock kernel launcher
macro cuda(ex)
# In real CUDA, this compiles kernel for GPU
# Type-unstable code would fail here!
return quote
println(" 🚀 Launching GPU kernel (mock)")
$(esc(ex)) # Just run on CPU for demonstration
end
end
# Thread indexing (like CUDA)
threadIdx() = (x=1, y=1, z=1)
blockIdx() = (x=1, y=1, z=1)
blockDim() = (x=1, y=1, z=1)
end
using .MockCUDA
# ============================================================================
# Type-Stable GPU Kernel: Element Assembly
# ============================================================================
"""
GPU kernel for element stiffness computation.
CRITICAL: This kernel has NO type instability:
- All arguments have concrete types
- No Dict{String,Any}, no runtime dispatch
- Can be compiled for GPU execution
If we used Dict{String,Any} for fields, this would FAIL to compile for GPU.
"""
function assemble_element_kernel!(
K_elements::CuArray{Float64,2}, # Pre-allocated output (n_elements, 64)
nodes::CuArray{Float64,2}, # Node coordinates (3, n_nodes)
connectivity::CuArray{Int,2}, # Element connectivity (8, n_elements)
E::Float64, # Young's modulus (type-stable!)
ν::Float64, # Poisson's ratio (type-stable!)
n_elements::Int
)
# GPU thread indexing (in real CUDA, this runs on GPU threads)
idx = MockCUDA.threadIdx().x +
(MockCUDA.blockIdx().x - 1) * MockCUDA.blockDim().x
if idx <= n_elements
# Extract element nodes (type-stable access)
elem_nodes = connectivity.data[:, idx]
# Mock stiffness computation (simplified)
# In reality, this would integrate over gauss points
K_local = E * (1 - ν^2) # Simplified scalar for demonstration
# Store result (in reality, this would be 8x8 matrix)
K_elements.data[idx, 1] = K_local
end
return nothing
end
# ============================================================================
# GPU Demonstration
# ============================================================================
println("Part 1: GPU Data Transfer and Kernel Execution")
println("-"^70)
# Setup problem data (type-stable!)
n_nodes = 1000
n_elements = 100
nodes = rand(Float64, 3, n_nodes) # Typed array: 3D coordinates
connectivity = rand(1:n_nodes, 8, n_elements) # Typed array: element topology
E = 210e3 # Concrete type: Float64
ν = 0.3 # Concrete type: Float64
println("\n✓ Created typed data structures:")
println(" - nodes: Array{Float64,2} ($(size(nodes)))")
println(" - connectivity: Array{Int,2} ($(size(connectivity)))")
println(" - E: Float64 = $E")
println(" - ν: Float64 = $ν")
# Transfer to GPU
println("\n✓ Transferring data to GPU:")
d_nodes = MockCUDA.cu(nodes)
d_connectivity = MockCUDA.cu(connectivity)
d_K_elements = MockCUDA.cu(zeros(Float64, n_elements, 64))
# Launch kernel
println("\n✓ Launching GPU kernel:")
MockCUDA.@cuda assemble_element_kernel!(
d_K_elements, d_nodes, d_connectivity, E, ν, n_elements
)
# Transfer results back
println("\n✓ Transferring results from GPU:")
K_elements = Array(d_K_elements)
println("\n✅ GPU execution successful!")
println(" Key insight: Type-stable data (Float64, Matrix{Float64}) transfers")
println(" directly to GPU with fast memcpy. No serialization needed.")
println()
# ============================================================================
# Mock MPI Interface (Minimal MPI-like interface without dependency)
# ============================================================================
"""
Mock MPI interface demonstrating efficient data transfer patterns.
In real MPI.jl:
MPI.Send(data, dest, tag, comm) # Uppercase = typed buffer transfer
MPI.send(data, dest, tag, comm) # Lowercase = slow serialization
"""
module MockMPI
struct Comm
rank::Int
size::Int
end
COMM_WORLD = Comm(0, 2)
function Send(data::Array{T,N}, dest::Int, tag::Int, comm::Comm) where {T,N}
nbytes = sizeof(data)
println(" 📨 MPI.Send: $(nbytes) bytes of $(eltype(data)) to rank $dest (fast buffer transfer)")
return nbytes
end
function Recv!(data::Array{T,N}, source::Int, tag::Int, comm::Comm) where {T,N}
nbytes = sizeof(data)
println(" 📬 MPI.Recv: $(nbytes) bytes of $(eltype(data)) from rank $source (fast buffer transfer)")
return nbytes
end
function send(data::Any, dest::Int, tag::Int, comm::Comm)
println(" 📨 MPI.send: serializing $(typeof(data)) to rank $dest (SLOW!)")
println(" ⚠️ Warning: This is ~100× slower than typed buffer transfer")
return 0
end
function recv(source::Int, tag::Int, comm::Comm)
println(" 📬 MPI.recv: deserializing from rank $source (SLOW!)")
return nothing
end
end
using .MockMPI
# ============================================================================
# MPI Demonstration
# ============================================================================
println("Part 2: MPI Data Transfer Patterns")
println("-"^70)
comm = MockMPI.COMM_WORLD
rank = comm.rank
size = comm.size
println("\n✓ MPI Communicator: rank=$rank, size=$size")
# Type-stable data transfer (FAST)
println("\n✓ Fast transfer: Typed arrays (uppercase MPI.Send)")
displacement = rand(Float64, 3, n_nodes)
forces = rand(Float64, 3, n_nodes)
MockMPI.Send(displacement, 1, 0, comm) # Uppercase = fast
MockMPI.Send(forces, 1, 1, comm)
# Type-unstable data transfer (SLOW)
println("\n✗ Slow transfer: Mixed types (lowercase MPI.send)")
fields_dict = Dict{String,Any}(
"displacement" => displacement,
"E" => E,
"nu" => ν
)
MockMPI.send(fields_dict, 1, 2, comm) # Lowercase = slow serialization
println("\n✅ MPI demonstration complete!")
println(" Key insight: Typed arrays (Matrix{Float64}) transfer ~100× faster")
println(" than mixed-type dictionaries (Dict{String,Any}).")
println()
# ============================================================================
# Combined GPU + MPI Pattern
# ============================================================================
println("Part 3: Combined GPU + MPI Workflow")
println("-"^70)
println("\n✓ Typical distributed GPU computation:")
println(" 1. Each MPI rank owns a subdomain")
println(" 2. Subdomain data (typed!) transfers to GPU")
println(" 3. GPU computes local contribution")
println(" 4. Results transfer back to CPU")
println(" 5. MPI exchanges boundary data (typed!)")
# Simulate subdomain on this rank
subdomain_nodes = rand(Float64, 3, n_nodes ÷ size)
subdomain_connectivity = rand(1:(n_nodes÷size), 8, n_elements ÷ size)
println("\n✓ Rank $rank subdomain:")
println(" - nodes: $(size(subdomain_nodes))")
println(" - elements: $(size(subdomain_connectivity, 2))")
# Transfer subdomain to GPU
println("\n✓ Transfer subdomain to GPU:")
d_sub_nodes = MockCUDA.cu(subdomain_nodes)
d_sub_connectivity = MockCUDA.cu(subdomain_connectivity)
d_sub_K = MockCUDA.cu(zeros(Float64, size(subdomain_connectivity, 2), 64))
# Compute on GPU
println("\n✓ Compute on GPU:")
MockCUDA.@cuda assemble_element_kernel!(
d_sub_K, d_sub_nodes, d_sub_connectivity, E, ν, size(subdomain_connectivity, 2)
)
# Transfer results back
println("\n✓ Transfer results from GPU:")
sub_K = Array(d_sub_K)
# Exchange boundary data with neighbor ranks
println("\n✓ MPI exchange boundary data:")
boundary_displacements = rand(Float64, 3, 10) # Mock boundary nodes
MockMPI.Send(boundary_displacements, (rank + 1) % size, 10, comm)
received_buffer = zeros(Float64, 3, 10)
MockMPI.Recv!(received_buffer, (rank - 1 + size) % size, 10, comm)
println("\n✅ Combined GPU+MPI workflow complete!")
println()
# ============================================================================
# Summary and Key Insights
# ============================================================================
println("="^70)
println("SUMMARY: Why Type Stability Matters for GPU/MPI")
println("="^70)
println("""
1. GPU Execution:
Type-stable code (Float64, Matrix{Float64}) compiles for GPU
Type-unstable code (Any, Dict{String,Any}) FAILS to compile
Transfer speed: ~1 GB/s for typed arrays (fast memcpy)
2. MPI Communication:
Typed arrays: MPI.Send (uppercase) = fast buffer transfer
Mixed types: MPI.send (lowercase) = slow serialization
Speed difference: ~100× faster for typed arrays
3. Zero Allocations:
Pre-allocated buffers on GPU/CPU
No allocations in kernel (required for GPU)
In-place operations preserve type stability
4. The Pattern:
- Define typed data structures (Matrix{Float64}, not Dict{String,Any})
- Pre-allocate buffers (cache, output arrays)
- Write type-stable kernels/functions
- Same code works on CPU, GPU, and across MPI
5. Performance Impact:
- CPU: 9-92× speedup (measured in field_storage_comparison.jl)
- GPU: Enables execution (type-unstable code cannot compile)
- MPI: 100× faster transfer (typed vs serialized)
CRITICAL INSIGHT:
Type stability is not a CPU optimization—it's a REQUIREMENT for GPU and
efficient MPI. The v0.5.1 Dict{String,Any} pattern makes GPU execution
impossible and MPI communication slow.
Any v1.0 design must ensure type-stable field access, regardless of where
data is stored (elements, global arrays, or elsewhere).
""")
println("="^70)
println()
-586
View File
@@ -1,586 +0,0 @@
#!/usr/bin/env julia
#
# GPU NodeSet MatVec Demo: NODAL ASSEMBLY for Krylov Methods
#
# Key insights:
# 1. JuliaFEM uses NODAL ASSEMBLY, not element assembly!
# 2. Loop over nodes, gather from connected elements
# 3. For Krylov (CG, GMRES), we need y = K*x (matvec)
# 4. Fields accessed through node_set (not passed separately)
# 5. This is TRULY general - no manual parameter extraction!
#
# Pattern:
# for node in node_set.nodes
# local_dofs = get_dofs(node)
# # Gather from all elements connected to this node
# y_local = compute_nodal_contribution(node, node_set.fields, x)
# y[local_dofs] = y_local # Direct write (no atomics!)
# end
using LinearAlgebra
using Printf
println("="^70)
println("GPU ElementSet MatVec Demo: Krylov-Ready Pattern")
println("="^70)
println()
# ============================================================================
# Mock Structures (NODAL ASSEMBLY)
# ============================================================================
"""Node with position"""
struct Node
id::UInt
x::Float64
y::Float64
z::Float64
end
"""Element references nodes (for gathering)"""
struct Element{N,B}
id::UInt
connectivity::NTuple{N,UInt} # Node IDs
basis::B
end
"""Mock basis type"""
struct MockBasis end
"""
NodeSet: Groups nodes + fields + connectivity
CRITICAL:
- Fields live at NODES (not elements!)
- node_to_elements[i] = list of elements connected to node i
- This enables nodal assembly: loop over nodes, gather from connected elements
"""
struct NodeSet{F}
name::String
nodes::Vector{Node}
elements::Vector{Element{4,MockBasis}} # For gathering
node_to_elements::Vector{Vector{Int}} # Inverse connectivity
fields::F # Type-stable field container (nodal fields!)
end
# Helper: Get DOF indices for a node (assuming 3 DOF per node)
function get_dofs(node::Node, dofs_per_node::Int=3)
return tuple(UInt.((node.id - 1) * dofs_per_node .+ (1:dofs_per_node))...)
end
# ============================================================================
# Mock GPU Module
# ============================================================================
module MockGPU
struct CuArray{T,N}
data::Array{T,N}
end
Base.length(a::CuArray) = length(a.data)
Base.getindex(a::CuArray, i...) = getindex(a.data, i...)
Base.setindex!(a::CuArray, v, i...) = setindex!(a.data, v, i...)
cu(x::Array) = CuArray(x)
cu(x::Vector) = CuArray(x)
cpu(x::CuArray) = x.data
macro cuda(args...)
func_call = args[end]
return esc(quote
$func_call
end)
end
export CuArray, cu, cpu, @cuda
end
using .MockGPU
# ============================================================================
# GPU Kernel: Matrix-Free Matrix-Vector Product (NODAL ASSEMBLY!)
# ============================================================================
"""
GPU kernel for matrix-vector product: y = K*x (NODAL ASSEMBLY)
TRULY GENERAL approach:
- Takes NodeSet (contains nodes + fields + connectivity)
- Loop over NODES (not elements!)
- Each node gathers from its connected elements
- No race conditions (each node writes to its own DOFs!)
- Accesses fields via node_set.fields (no manual extraction!)
- Computes nodal contribution by gathering from all connected elements
This is what Krylov methods need - NOT the full K matrix!
"""
function gpu_matvec_kernel!(
y::CuArray{Float64,1}, # Output: y = K*x
x::CuArray{Float64,1}, # Input vector
node_set::NodeSet, # Contains nodes + fields + connectivity!
dofs_per_node::Int,
)
# In real CUDA: thread_id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
# Each thread processes one NODE
for node_id in 1:length(node_set.nodes)
node = node_set.nodes[node_id]
# Access fields from node_set (GENERAL!)
E = node_set.fields.E
ν = node_set.fields.ν
# Get this node's DOFs
local_dofs = get_dofs(node, dofs_per_node)
n_local_dofs = length(local_dofs)
# Initialize nodal contribution to zero
y_nodal = ntuple(i -> 0.0, n_local_dofs)
# Gather from all elements connected to this node
connected_elements = node_set.node_to_elements[node_id]
for elem_idx in connected_elements
element = node_set.elements[elem_idx]
# Find this node's position in element connectivity
local_node_idx = findfirst(==(node.id), element.connectivity)
# Get all element DOFs
elem_dofs = UInt[]
for node_id_in_elem in element.connectivity
for d in 1:dofs_per_node
push!(elem_dofs, (node_id_in_elem - 1) * dofs_per_node + d)
end
end
# Extract element x values
x_elem = [x[dof] for dof in elem_dofs]
# Mock element stiffness contribution
# In real code: K_elem = assemble_element_stiffness(element, E, ν)
# Extract rows corresponding to this node
K_elem_factor = E * (1 - ν^2) * 0.1
# Add this element's contribution to nodal y
# (rows corresponding to this node)
node_start = (local_node_idx - 1) * dofs_per_node
for i in 1:n_local_dofs
row_in_elem = node_start + i
# Sum over all DOFs in element
for j in 1:length(x_elem)
y_nodal = ntuple(k -> k == i ? y_nodal[k] + K_elem_factor * x_elem[j] : y_nodal[k], n_local_dofs)
end
end
end
# Write nodal contribution to global y (no atomics needed!)
for i in 1:n_local_dofs
y[local_dofs[i]] = y_nodal[i]
end
end
return nothing
end
# ============================================================================
# CPU Version: Same Logic
# ============================================================================
"""CPU matrix-vector product using nodal assembly pattern"""
function cpu_matvec!(
y::Vector{Float64},
x::Vector{Float64},
node_set::NodeSet,
dofs_per_node::Int,
)
# Zero output
fill!(y, 0.0)
# Access fields from node_set (GENERAL!)
E = node_set.fields.E
ν = node_set.fields.ν
# Loop over NODES (not elements!)
for (node_id, node) in enumerate(node_set.nodes)
# Get this node's DOFs
local_dofs = get_dofs(node, dofs_per_node)
n_local_dofs = length(local_dofs)
# Initialize nodal contribution
y_nodal = zeros(n_local_dofs)
# Gather from all connected elements
connected_elements = node_set.node_to_elements[node_id]
for elem_idx in connected_elements
element = node_set.elements[elem_idx]
# Find this node's position in element
local_node_idx = findfirst(==(node.id), element.connectivity)
# Get all element DOFs
elem_dofs = Int[]
for node_id_in_elem in element.connectivity
for d in 1:dofs_per_node
push!(elem_dofs, (node_id_in_elem - 1) * dofs_per_node + d)
end
end
# Extract element x values
x_elem = [x[dof] for dof in elem_dofs]
# Mock element stiffness computation
K_elem_factor = E * (1 - ν^2) * 0.1
# Add element contribution (rows for this node)
node_start = (local_node_idx - 1) * dofs_per_node
for i in 1:n_local_dofs
for j in 1:length(x_elem)
y_nodal[i] += K_elem_factor * x_elem[j]
end
end
end
# Write nodal contribution to global y
for i in 1:n_local_dofs
y[local_dofs[i]] = y_nodal[i]
end
end
return y
end
# ============================================================================
# Setup Problem (NODAL ASSEMBLY)
# ============================================================================
println("Setting up problem with NODAL ASSEMBLY...")
println()
# Create simple 2D mesh
n_x = 10
n_y = 10
n_nodes = n_x * n_y
dofs_per_node = 3
n_dofs = n_nodes * dofs_per_node
# Create nodes
nodes = [
Node(
UInt((j-1)*n_x + i),
Float64(i),
Float64(j),
0.0
)
for j in 1:n_y for i in 1:n_x
]
# Create elements (quads)
elements = Element{4,MockBasis}[]
for j in 1:(n_y-1)
for i in 1:(n_x-1)
node1 = UInt((j-1)*n_x + i)
node2 = UInt((j-1)*n_x + i + 1)
node3 = UInt(j*n_x + i + 1)
node4 = UInt(j*n_x + i)
push!(elements, Element{4,MockBasis}(
UInt(length(elements) + 1),
(node1, node2, node3, node4),
MockBasis()
))
end
end
n_elements = length(elements)
# Build node-to-elements connectivity
node_to_elements = [Int[] for _ in 1:n_nodes]
for (elem_idx, element) in enumerate(elements)
for node_id in element.connectivity
push!(node_to_elements[node_id], elem_idx)
end
end
# Fields in node set (type-stable!)
fields = (
E=210e3,
ν=0.3,
)
# Create node set (THIS IS THE REAL STRUCTURE!)
node_set = NodeSet("steel_body", nodes, elements, node_to_elements, fields)
println("Problem setup:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println(" DOFs per node: $dofs_per_node")
println(" Total DOFs: $n_dofs")
println(" Node type: ", typeof(nodes[1]))
println(" Field type: ", typeof(node_set.fields))
println(" Fields accessed via: node_set.fields.E, node_set.fields.ν")
println(" Average elements/node: ", sum(length.(node_to_elements)) / n_nodes)
println()
println("NODAL ASSEMBLY:")
println(" - Loop over nodes (not elements!)")
println(" - Each node gathers from connected elements")
println(" - No race conditions (each node owns its DOFs)")
println()
# ============================================================================
# Test CPU MatVec
# ============================================================================
println("="^70)
println("CPU Matrix-Vector Product")
println("="^70)
println()
x = randn(n_dofs)
y_cpu = zeros(n_dofs)
println("Computing y = K*x using NODAL ASSEMBLY pattern...")
@time cpu_matvec!(y_cpu, x, node_set, dofs_per_node)
println("\nResult:")
println(" ||x||: ", @sprintf("%.6e", norm(x)))
println(" ||y||: ", @sprintf("%.6e", norm(y_cpu)))
println(" ✓ Matrix-vector product complete (nodal assembly)")
println()
# ============================================================================
# Test GPU MatVec (Mock)
# ============================================================================
println("="^70)
println("GPU Matrix-Vector Product (Mock CUDA) - NODAL ASSEMBLY")
println("="^70)
println()
println("Key insight: NodeSet goes to GPU!")
println(" - Nodes: node_set.nodes")
println(" - Elements: node_set.elements (for gathering)")
println(" - Connectivity: node_set.node_to_elements")
println(" - Fields: node_set.fields")
println(" - No manual parameter extraction needed!")
println(" - Each GPU thread processes ONE NODE (not element)")
println()
# Transfer to GPU
x_gpu = cu(x)
y_gpu = cu(zeros(n_dofs))
println("Launching GPU kernel (one thread per NODE)...")
@cuda threads = 256 blocks = ceil(Int, n_nodes / 256) gpu_matvec_kernel!(
y_gpu, x_gpu, node_set, dofs_per_node
)
println(" ✓ Kernel execution complete")
println()
# Transfer back and verify
y_gpu_result = cpu(y_gpu)
error = norm(y_gpu_result - y_cpu) / (norm(y_cpu) + 1e-10)
println("📊 GPU Results:")
println(" ||y_GPU||: ", @sprintf("%.6e", norm(y_gpu_result)))
println(" Relative error: ", @sprintf("%.6e", error))
println(" ✓ GPU matches CPU: ", error < 1e-6 ? "YES ✅" : "NO ❌")
println()
println("🎯 NODAL ASSEMBLY ADVANTAGES:")
println(" ✅ No atomic operations needed (each node owns its DOFs)")
println(" ✅ Natural for contact mechanics (contact forces at nodes)")
println(" ✅ Clean domain decomposition (node ownership)")
println(" ✅ Better cache locality (node data grouped)")
println()
# ============================================================================
# Demonstrate GMRES Pattern (Conceptual)
# ============================================================================
println("="^70)
println("How This Enables Krylov Methods (GMRES/CG)")
println("="^70)
println()
println("""
GMRES only needs matrix-vector products, not the matrix itself!
Traditional approach (WRONG for large problems):
K = assemble_global_matrix(elements) # O(N²) memory!
y = K * x # Dense operation
Matrix-free approach (CORRECT):
y = matvec(element_set, x) # O(N) memory!
# Computed by looping over elements, no global K
GMRES iteration:
for iteration in 1:max_iterations
# Build Krylov subspace using matvec
v_new = matvec(element_set, v_old) # ← Our GPU kernel!
# Orthogonalize (Arnoldi)
# ...
# Check convergence
if residual < tolerance
break
end
end
For contact mechanics + plasticity:
- Update contact state (node-by-node)
- Update material state (integration points)
- Compute y = K_tangent * x using current state
- No need to form K_tangent explicitly!
Our GPU kernel is PERFECT for this:
1. Element-local computation (easy to parallelize)
2. Fields accessed naturally (element_set.fields)
3. Returns y vector (what GMRES needs)
4. O(N) memory (no global matrix)
""")
# ============================================================================
# Show Field Access Pattern
# ============================================================================
println("="^70)
println("Field Access Pattern (TRULY General)")
println("="^70)
println()
println("""
In the GPU kernel, we access fields like this:
function gpu_matvec_kernel!(y, x, element_set, ...)
element = element_set.elements[elem_id]
# Access fields from element_set (not passed separately!)
E = element_set.fields.E # ← GENERAL!
ν = element_set.fields.ν # ← GENERAL!
u = element_set.fields.u # ← If displacement field exists
# Get local DOFs from element connectivity
local_dofs = get_dofs(element)
# Extract local portion of x
x_local = x[local_dofs]
# Compute local matvec
y_local = K_local(element, E, ν) * x_local
# Add to global (atomic on GPU)
y[local_dofs] += y_local
end
No manual parameter extraction!
No separate arrays for E, ν, etc!
Everything accessed through element_set!
This works for ANY field type (NamedTuple, struct, whatever) as long as
it's type-stable!
""")
# ============================================================================
# Time-Dependent Fields Example
# ============================================================================
println("="^70)
println("Time-Dependent Fields (Bonus)")
println("="^70)
println()
println("For transient problems, create new node_set each time step:")
println()
println("""
# Time stepping loop
for t in time_steps
# Solve for new displacement using GMRES (matrix-free!)
u_new = gmres(x0) do x
y = zeros(n_dofs)
gpu_matvec_kernel!(y, x, node_set, dofs_per_node)
return y
end
# Create NEW field container (cheap!)
fields_new = (
E = node_set.fields.E, # Constant (keep)
ν = node_set.fields.ν, # Constant (keep)
u = u_new, # Updated!
temperature = T_new, # Updated!
)
# Create new node set (cheap - just wraps references)
node_set = NodeSet(name, nodes, elements, node_to_elements, fields_new)
# Next iteration uses updated fields automatically!
end
Creating new NamedTuple: ~2-3 ns (just wraps references)
No data copying needed!
""")
# ============================================================================
# Summary
# ============================================================================
println("="^70)
println("SUMMARY: Matrix-Free Krylov with NODAL ASSEMBLY")
println("="^70)
println()
println("""
TRULY GENERAL approach with NODAL ASSEMBLY:
1. NodeSet contains nodes + elements + connectivity + fields
2. Fields accessed via node_set.fields (no manual extraction!)
3. GPU kernel loops over NODES (not elements!)
4. Each node gathers from connected elements
5. No race conditions (each node owns its DOFs)
6. O(N) memory (no global matrix)
For GMRES/CG:
- Only need matvec operation (y = K*x)
- No need to form or store K
- Perfect for contact + plasticity (state-dependent K)
- Scales to millions of DOFs
Nodal assembly advantages:
- node_set.fields.E (constant material property)
- node_set.fields.ν (constant material property)
- node_set.fields.u (nodal displacement)
- Contact forces natural (at nodes!)
- Domain decomposition clean (node ownership)
- No atomic operations needed on GPU
Type stability:
- NodeSet{F} has known field type F
- Compiler generates optimal code
- Zero runtime dispatch
- GPU-compatible
This is the NODAL ASSEMBLY pattern for JuliaFEM v1.0!
🎯 Why nodal assembly?
1. Contact mechanics is nodal (forces, constraints at nodes)
2. Domain decomposition is nodal (clean node ownership)
3. No atomic operations on GPU (each node owns its DOFs)
4. Better cache locality (node data grouped together)
5. Natural for adaptive refinement (local node operations)
""")
println("="^70)
println("✅ DEMONSTRATION COMPLETE - NODAL ASSEMBLY")
println("="^70)
println()
println("Next steps:")
println(" 1. Implement real nodal matvec with full stiffness computation")
println(" 2. Integrate with Krylov.jl for GMRES/CG")
println(" 3. Add contact state updates (natural at nodes!)")
println(" 4. Add material state updates (at integration points)")
println(" 5. Test on real CUDA hardware")
println(" 6. Benchmark: nodal vs element assembly")
println()
-203
View File
@@ -1,203 +0,0 @@
#!/usr/bin/env julia
#
# GPU-Only Demonstration (No MPI Required)
#
# Demonstrates that type-stable field data can execute on real CUDA GPU.
# Simpler than full GPU+MPI demo - just shows GPU capability.
#
using LinearAlgebra
println("="^70)
println("GPU Type-Stable Kernel Demonstration")
println("="^70)
println()
# Try to load CUDA
CUDA_AVAILABLE = false
try
@eval using CUDA
if CUDA.functional()
global CUDA_AVAILABLE = true
println("✅ CUDA GPU detected: $(CUDA.name(CUDA.device()))")
println(" Memory: $(CUDA.totalmem(CUDA.device()) ÷ 10^9) GB")
println()
else
println("❌ CUDA.jl loaded but no functional GPU detected")
exit(1)
end
catch e
println("❌ CUDA.jl not available: $e")
println(" Install with: using Pkg; Pkg.add(\"CUDA\")")
exit(1)
end
println("="^70)
println("Demonstration: Type-Stable Assembly Kernel")
println("="^70)
println()
# Problem setup
n_nodes = 10000
n_elements = 1000
println("Setup:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elements")
println()
# Type-stable data structures
nodes = rand(Float64, 3, n_nodes)
connectivity = rand(1:n_nodes, 8, n_elements)
E = 210e3 # Young's modulus
ν = 0.3 # Poisson's ratio
println("Data types (type-stable):")
println(" nodes: $(typeof(nodes))")
println(" connectivity: $(typeof(connectivity))")
println(" E: $(typeof(E))")
println(" ν: $(typeof(ν))")
println()
# Define GPU kernel
function assemble_element_kernel!(
K_elements::CuDeviceMatrix{Float64},
nodes::CuDeviceMatrix{Float64},
connectivity::CuDeviceMatrix{Int32},
E::Float64,
ν::Float64,
n_elements::Int32
)
# GPU thread indexing
idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if idx <= n_elements
# Type-stable access to element data
# In real FEM: would integrate over Gauss points
# Here: simplified computation to demonstrate GPU execution
# Mock stiffness calculation
K_local = E * (1 - ν^2)
# Store result
K_elements[idx, 1] = K_local
end
return nothing
end
println("="^70)
println("Step 1: Transfer Data to GPU")
println("="^70)
println()
# Calculate sizes
nodes_bytes = sizeof(nodes)
conn_bytes = sizeof(connectivity)
total_bytes = nodes_bytes + conn_bytes
println("Transferring to GPU:")
println(" nodes: $(nodes_bytes ÷ 1024) KB")
println(" connectivity: $(conn_bytes ÷ 1024) KB")
println(" Total: $(total_bytes ÷ 1024) KB")
println()
# Transfer to GPU
d_nodes = CuArray(nodes)
d_connectivity = CuArray(Int32.(connectivity))
d_K_elements = CUDA.zeros(Float64, n_elements, 64)
println("✅ Data on GPU")
println()
println("="^70)
println("Step 2: Launch GPU Kernel")
println("="^70)
println()
# Kernel launch configuration
threads_per_block = 256
blocks = cld(n_elements, threads_per_block)
println("Kernel configuration:")
println(" Threads per block: $threads_per_block")
println(" Blocks: $blocks")
println(" Total threads: $(blocks * threads_per_block)")
println()
println("Launching kernel...")
@cuda threads = threads_per_block blocks = blocks assemble_element_kernel!(
d_K_elements, d_nodes, d_connectivity, E, ν, Int32(n_elements)
)
# Wait for completion
CUDA.synchronize()
println("✅ Kernel executed successfully")
println()
println("="^70)
println("Step 3: Transfer Results from GPU")
println("="^70)
println()
K_elements = Array(d_K_elements)
result_bytes = sizeof(K_elements)
println("Transferred from GPU:")
println(" Results: $(result_bytes ÷ 1024) KB")
println()
println("="^70)
println("Step 4: Verify Results")
println("="^70)
println()
expected_value = E * (1 - ν^2)
actual_values = K_elements[:, 1]
all_match = all(abs.(actual_values .- expected_value) .< 1e-10)
println("Verification:")
println(" Expected value: $expected_value")
println(" First result: $(actual_values[1])")
println(" All elements match: $(all_match ? "" : "")")
println()
if all_match
println("="^70)
println("✅ SUCCESS: Type-Stable GPU Kernel Executed Correctly")
println("="^70)
println()
println("Key Achievements:")
println()
println("1. Type-stable kernel compiled for GPU")
println(" - All parameters have concrete types (Float64, Int32)")
println(" - No Dict{String,Any} or runtime dispatch")
println(" - Compiler generated optimized GPU machine code")
println()
println("2. Fast GPU memory transfer")
println(" - Typed arrays transferred as contiguous buffers")
println(" - No serialization overhead")
println(" - Same pattern works for MPI communication")
println()
println("3. Zero allocations in kernel")
println(" - All arrays pre-allocated")
println(" - In-place operations only")
println(" - Required for GPU execution")
println()
println("Why This Matters:")
println()
println("• Dict{String,Any} field storage CANNOT compile for GPU")
println(" - Compiler error: cannot determine types")
println(" - Would prevent any GPU acceleration")
println()
println("• Type-stable storage (Matrix{Float64}) works everywhere:")
println(" - CPU: 9-92× faster (measured)")
println(" - GPU: Enables execution (demonstrated)")
println(" - MPI: Fast transfers (same pattern)")
println()
println("Conclusion: Type stability is not optional for modern HPC")
println("="^70)
else
println("❌ FAILED: Results don't match expected values")
end
-400
View File
@@ -1,400 +0,0 @@
#!/usr/bin/env julia
#
# Multi-GPU MPI Krylov Solver Demonstration
#
# This script demonstrates a distributed FEM-like solver using:
# 1. Nodal assembly pattern (row-by-row matrix construction)
# 2. Multi-GPU setup with MPI communication
# 3. Krylov iterative solver (Conjugate Gradient)
# 4. Verification against exact solution
#
# Requirements:
# - CUDA-capable GPU on each MPI rank (optional, will use CPU if unavailable)
# - MPI installation
# - Run with: mpiexec -np 2 julia --project=. benchmarks/krylov_mpi_gpu_demo.jl
#
# KEY INSIGHT: Nodal assembly + type stability enables distributed GPU solving
#
using LinearAlgebra
using Random
using Printf
# Try to load CUDA (optional, will fall back to CPU)
CUDA_AVAILABLE = false
try
using CUDA
if CUDA.functional()
global CUDA_AVAILABLE = true
println("✓ CUDA GPU detected on this rank: $(CUDA.name(CUDA.device()))")
else
println("⚠ CUDA.jl installed but no GPU detected on this rank")
end
catch e
println(" CUDA.jl not available (will use CPU): $e")
end
# Load MPI (required)
using MPI
MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
nranks = MPI.Comm_size(comm)
# Helper for master-only printing
function println_master(args...)
if rank == 0
println(args...)
end
end
println_master("="^70)
println_master("Multi-GPU MPI Krylov Solver Demonstration")
println_master("="^70)
println_master("Configuration:")
println_master(" MPI ranks: $nranks")
println_master(" CUDA available: $CUDA_AVAILABLE")
println_master("")
#==============================================================================
Part 1: Generate Test Problem (Distributed Nodal Assembly)
==============================================================================#
println_master("Part 1: Generating Test Problem")
println_master("-"^70)
# Problem size (total DOFs)
const N = 10
println_master(" Problem size: $(N)×$(N) system")
# Each rank owns a partition of nodes (rows)
nodes_per_rank = div(N, nranks)
remainder = N % nranks
my_start = rank * nodes_per_rank + min(rank, remainder) + 1
my_end = my_start + nodes_per_rank - 1 + (rank < remainder ? 1 : 0)
my_n = my_end - my_start + 1
println(" Rank $rank: owns nodes $(my_start):$(my_end) ($(my_n) nodes)")
# Generate global problem (same on all ranks for verification)
# In real FEM: each rank would only know local + ghost nodes
Random.seed!(12345) # Same seed on all ranks for reproducibility
# Create a symmetric positive definite matrix (fake assembly)
# In real FEM: this would come from element integration
A_global = rand(Float64, N, N)
A_global = A_global' * A_global # Make SPD
A_global += 10.0 * I(N) # Ensure strong diagonal dominance
# Exact solution (known)
x_exact = Float64[i for i in 1:N]
# Right-hand side
b_global = A_global * x_exact
println_master(" ✓ Generated SPD matrix (condition number ≈ $(cond(A_global)))")
println_master(" ✓ Exact solution: x = [1, 2, 3, ..., $N]")
println_master("")
#==============================================================================
Part 2: Nodal Assembly Pattern (Row-by-Row)
==============================================================================#
println_master("Part 2: Nodal Assembly Pattern")
println_master("-"^70)
"""
get_row(A_global, i::Int) -> Vector{Float64}
Simulate nodal assembly: returns the i-th row of the system matrix.
In real FEM: this would assemble contributions from all elements
connected to node i.
"""
function get_row(A_global, i::Int)
return A_global[i, :]
end
"""
get_rhs(b_global, i::Int) -> Float64
Get right-hand side value for node i.
"""
function get_rhs(b_global, i::Int)
return b_global[i]
end
# Each rank assembles its local rows
my_rows = Matrix{Float64}(undef, my_n, N)
my_rhs = Vector{Float64}(undef, my_n)
for (local_i, global_i) in enumerate(my_start:my_end)
my_rows[local_i, :] = get_row(A_global, global_i)
my_rhs[local_i] = get_rhs(b_global, global_i)
end
println(" Rank $rank: assembled $(my_n) rows locally")
println_master(" ✓ Nodal assembly complete (each rank has its partition)")
println_master("")
#==============================================================================
Part 3: GPU Transfer (Optional, if CUDA available)
==============================================================================#
if CUDA_AVAILABLE
println_master("Part 3: GPU Transfer")
println_master("-"^70)
# Transfer local data to GPU
d_my_rows = CuArray(my_rows)
d_my_rhs = CuArray(my_rhs)
bytes_transferred = sizeof(my_rows) + sizeof(my_rhs)
println(" Rank $rank: transferred $(bytes_transferred) bytes to GPU")
println_master(" ✓ Each rank transferred local data to its GPU")
println_master("")
else
println_master("Part 3: GPU Transfer")
println_master("-"^70)
println_master(" ⚠ CUDA not available, using CPU arrays")
println_master("")
d_my_rows = my_rows
d_my_rhs = my_rhs
end
#==============================================================================
Part 4: Distributed Matrix-Vector Product
==============================================================================#
println_master("Part 4: Distributed Matrix-Vector Product")
println_master("-"^70)
"""
matvec_distributed!(y_local, x_global, A_local)
Compute y_local = A_local * x_global (distributed matrix-vector product).
Each rank computes its portion of the result using its local rows.
"""
function matvec_distributed!(y_local, x_global, A_local)
# Copy to device if using GPU
if CUDA_AVAILABLE
d_x = CuArray(x_global)
d_A = isa(A_local, CuArray) ? A_local : CuArray(A_local)
d_y = d_A * d_x
copyto!(y_local, Array(d_y))
else
# CPU computation
mul!(y_local, A_local, x_global)
end
return nothing
end
# Test matvec
x_test = ones(Float64, N)
y_test = Vector{Float64}(undef, my_n)
matvec_distributed!(y_test, x_test, CUDA_AVAILABLE ? Array(d_my_rows) : my_rows)
println(" Rank $rank: matvec test complete ($(length(y_test)) outputs)")
println_master(" ✓ Distributed matrix-vector product working")
println_master("")
#==============================================================================
Part 5: Conjugate Gradient Solver (Distributed)
==============================================================================#
println_master("Part 5: Conjugate Gradient Solver")
println_master("-"^70)
"""
cg_distributed(A_local, b_local, x0; maxiter=100, tol=1e-6)
Distributed Conjugate Gradient solver.
Each rank holds local rows of A and local entries of vectors.
Uses MPI collectives for global dot products and norms.
"""
function cg_distributed(A_local, b_local, x0; maxiter=100, tol=1e-6)
n_global = length(x0)
n_local = length(b_local)
# Initial guess
x = copy(x0)
# Initial residual: r = b - A*x (distributed)
r_local = similar(b_local)
matvec_distributed!(r_local, x, A_local)
r_local .= b_local .- r_local
# Global residual norm
r_norm_local = dot(r_local, r_local)
r_norm_sq = MPI.Allreduce(r_norm_local, MPI.SUM, comm)
r_norm_0 = sqrt(r_norm_sq)
if rank == 0
@printf(" Initial residual: %.6e\n", r_norm_0)
end
# CG iteration
p = copy(x) # Search direction (global vector)
# Distribute initial r to all ranks for p initialization
# Each rank needs full vector for matvec
r_global = Vector{Float64}(undef, n_global)
# Gather r from all ranks
recvcounts = Int32.(MPI.Allgather(n_local, comm))
displs = Int32.([0; cumsum(recvcounts[1:end-1])])
MPI.Allgatherv!(r_local, r_global, recvcounts, comm)
p .= r_global
for iter in 1:maxiter
# Compute A*p (distributed)
Ap_local = similar(b_local)
matvec_distributed!(Ap_local, p, A_local)
# Global dot products: alpha = (r'*r) / (p'*A*p)
pAp_local = dot(r_local, r_local) # We stored r'*r from previous iteration
pAp_numerator = MPI.Allreduce(pAp_local, MPI.SUM, comm)
# Need p'*Ap - but p is global and Ap is local
# Gather Ap
Ap_global = Vector{Float64}(undef, n_global)
MPI.Allgatherv!(Ap_local, Ap_global, recvcounts, comm)
pAp_denominator = dot(p, Ap_global)
alpha = pAp_numerator / pAp_denominator
# Update solution and residual (global vectors)
x .+= alpha .* p
# Update local residual
r_local .-= alpha .* Ap_local
# Check convergence
r_norm_local = dot(r_local, r_local)
r_norm_sq = MPI.Allreduce(r_norm_local, MPI.SUM, comm)
r_norm = sqrt(r_norm_sq)
if rank == 0
@printf(" Iteration %3d: residual = %.6e (reduction: %.2f%%)\n",
iter, r_norm, 100.0 * (1.0 - r_norm / r_norm_0))
end
if r_norm < tol
if rank == 0
println(" ✓ Converged in $iter iterations")
end
return x, iter, r_norm
end
# Update search direction: beta = r_new'*r_new / r_old'*r_old
beta = r_norm_sq / pAp_numerator
# Gather updated r for next p
MPI.Allgatherv!(r_local, r_global, recvcounts, comm) # Update search direction
p .= r_global .+ beta .* p
end
if rank == 0
println(" ⚠ Did not converge in $maxiter iterations")
end
return x, maxiter, sqrt(r_norm_sq)
end
# Solve the system
x0 = zeros(Float64, N)
x_solution, iters, final_residual = cg_distributed(
CUDA_AVAILABLE ? Array(d_my_rows) : my_rows,
CUDA_AVAILABLE ? Array(d_my_rhs) : my_rhs,
x0,
maxiter=100,
tol=1e-10
)
println_master("")
#==============================================================================
Part 6: Verification Against Exact Solution
==============================================================================#
println_master("Part 6: Verification")
println_master("-"^70)
# Compute error
error = norm(x_solution - x_exact) / norm(x_exact)
println_master("Solution comparison:")
println_master(" Exact: ", join([@sprintf("%.3f", x) for x in x_exact[1:min(5, N)]], ", "),
N > 5 ? ", ..." : "")
println_master(" Computed: ", join([@sprintf("%.3f", x) for x in x_solution[1:min(5, N)]], ", "),
N > 5 ? ", ..." : "")
println_master("")
println_master(" Relative error: ", @sprintf("%.6e", error))
println_master(" Converged in: $iters iterations")
println_master(" Final residual: ", @sprintf("%.6e", final_residual))
println_master("")
if error < 1e-6
println_master("✅ VERIFICATION PASSED (error < 1e-6)")
else
println_master("❌ VERIFICATION FAILED (error = $error)")
end
println_master("")
#==============================================================================
Summary: Key Insights
==============================================================================#
println_master("="^70)
println_master("SUMMARY: Multi-GPU MPI Krylov Solver")
println_master("="^70)
println_master("")
println_master("✅ Demonstrated on Real Hardware:")
println_master("")
println_master("1. Nodal Assembly Pattern:")
println_master(" • Each rank assembles its local rows (nodes $(my_start):$(my_end))")
println_master(" • Row-by-row construction: get_row() abstraction")
println_master(" • Natural for contact mechanics (nodal basis)")
println_master("")
println_master("2. Distributed Computing:")
println_master(" • Problem split across $nranks MPI ranks")
println_master(" • Each rank owns $(nodes_per_rank) nodes")
println_master(" • MPI collectives for global operations (dot products)")
println_master("")
if CUDA_AVAILABLE
println_master("3. Multi-GPU Execution:")
println_master(" • Each rank transferred data to its local GPU")
println_master(" • Matrix-vector products computed on GPU")
println_master(" • Results synchronized via MPI")
println_master("")
else
println_master("3. CPU Execution:")
println_master(" • CUDA not available, used CPU arrays")
println_master(" • Same algorithm works on CPU/GPU")
println_master(" • Type stability enables both paths")
println_master("")
end
println_master("4. Krylov Iterative Solver:")
println_master(" • Conjugate Gradient (CG) method")
println_master(" • Distributed matrix-vector products")
println_master(" • Converged in $iters iterations")
println_master(" • Relative error: ", @sprintf("%.6e", error))
println_master("")
println_master("Key Insight:")
println_master(" Type-stable nodal assembly + distributed matvec → scalable solving")
println_master(" Same code pattern: CPU → GPU → MPI → Multi-GPU")
println_master("")
println_master("Relevance to JuliaFEM:")
println_master(" • Nodal assembly aligns with contact mechanics")
println_master(" • Row-by-row construction enables streaming assembly")
println_master(" • Type stability requirement validated on real hardware")
println_master(" • Distributed solving demonstrated at small scale")
println_master("")
println_master("="^70)
MPI.Finalize()
-513
View File
@@ -1,513 +0,0 @@
"""
Complete Newton-Krylov-Anderson Reference Implementation (CPU)
==============================================================
Shows the full nonlinear solver pipeline:
1. Newton iteration (outer loop)
2. GMRES for linear solve (inner loop, matrix-free)
3. Anderson acceleration (outer loop acceleration)
4. Perfect plasticity with state variables
This is the REFERENCE. Once it works, we port to GPU.
"""
using LinearAlgebra
using Tensors
using Printf
# ============================================================================
# Material: Perfect Plasticity with State Variables
# ============================================================================
struct VonMisesPlasticity
E::Float64 # Young's modulus
ν::Float64 # Poisson's ratio
σ_y::Float64 # Yield stress
end
# Material constants
λ(mat::VonMisesPlasticity) = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
μ(mat::VonMisesPlasticity) = mat.E / (2(1 + mat.ν))
# Plastic state per integration point
mutable struct PlasticState
ε_p::SymmetricTensor{2,3,Float64,6} # Plastic strain
α::Float64 # Accumulated plastic strain
end
PlasticState() = PlasticState(zero(SymmetricTensor{2,3,Float64}), 0.0)
"""
Compute stress with return mapping (radial return algorithm).
Returns (σ, state_new, plastic_loading).
"""
function compute_stress_with_plasticity(material::VonMisesPlasticity,
ε_total::SymmetricTensor{2,3,Float64},
state_old::PlasticState)
# Trial elastic strain
ε_e_trial = ε_total - state_old.ε_p
# Trial stress (elastic predictor)
λ_val = λ(material)
μ_val = μ(material)
I = one(ε_e_trial)
σ_trial = λ_val * tr(ε_e_trial) * I + 2 * μ_val * ε_e_trial
# Deviatoric stress
σ_dev = dev(σ_trial)
σ_eq = sqrt(3 / 2 * (σ_dev σ_dev)) # von Mises equivalent stress
# Check yield condition
f_trial = σ_eq - material.σ_y
if f_trial <= 0.0
# Elastic loading - no plasticity
return (σ_trial, state_old, false)
else
# Plastic loading - return mapping
Δγ = f_trial / (3 * μ_val) # Plastic multiplier (for perfect plasticity)
# Return mapping
n = σ_dev / σ_eq # Flow direction (normal to yield surface)
σ_new = σ_trial - 2 * μ_val * Δγ * n
# Update plastic strain
Δε_p = Δγ * n
ε_p_new = state_old.ε_p + Δε_p
α_new = state_old.α + Δγ
state_new = PlasticState(ε_p_new, α_new)
return (σ_new, state_new, true)
end
end
# ============================================================================
# Tet4: Linear Tetrahedron (simpler than Tet10 for GPU proof-of-concept)
# ============================================================================
# Gauss quadrature: 1-point for Tet4 (sufficient for linear element)
const GAUSS_TET4_1PT = (
(Vec{3}((0.25, 0.25, 0.25)), 1.0 / 6.0), # Weight = volume of reference tet
)
"""
Shape function derivatives for Tet4 (constant in reference element).
Returns tuple of 4 Vec{3}.
"""
function tet4_shape_derivatives()
# For reference Tet4: vertices at (0,0,0), (1,0,0), (0,1,0), (0,0,1)
# dN/dξ are constant:
dN1 = Vec{3}((-1.0, -1.0, -1.0))
dN2 = Vec{3}((1.0, 0.0, 0.0))
dN3 = Vec{3}((0.0, 1.0, 0.0))
dN4 = Vec{3}((0.0, 0.0, 1.0))
return (dN1, dN2, dN3, dN4)
end
# ============================================================================
# Element Residual Assembly (Matrix-Free)
# ============================================================================
"""
Assemble element residual for Tet4.
This is what gets called inside the matrix-free operator.
"""
function assemble_element_residual!(r_elem, X, u, material, states)
# Shape derivatives (constant for Tet4)
dN_dxi = tet4_shape_derivatives()
# Zero residual
fill!(r_elem, zero(Vec{3,Float64}))
# Integration loop (1 point for Tet4)
for (gp_idx, (xivec, w)) in enumerate(GAUSS_TET4_1PT)
# Jacobian: J = Σ dN_i ⊗ X_i
J = sum(dN_dxi[i] X[i] for i in 1:4)
detJ = det(J)
invJ = inv(J)
# Physical derivatives: dN/dx = J^-1 · dN/dξ
dN_dx = ntuple(i -> invJ Vec{3}(dN_dxi[i]), Val(4))
# Strain: ε = sym(∇u)
eps = symmetric(sum(dN_dx[i] u[i] for i in 1:4))
# Stress with plasticity
state_old = states[gp_idx]
(sigma, state_new, plastic) = compute_stress_with_plasticity(material, eps, state_old)
states[gp_idx] = state_new # Update state
# Nodal forces: f_i = dN_i · σ
f_contrib = ntuple(i -> dN_dx[i] sigma, Val(4))
# Accumulate
for i in 1:4
r_elem[i] += f_contrib[i] * (w * detJ)
end
end
return r_elem
end
"""
Global residual assembly (loop over elements).
"""
function assemble_global_residual!(r, u, mesh, material, element_states)
fill!(r, 0.0)
for (elem_idx, elem) in enumerate(mesh.elements)
# Extract element data
X = ntuple(i -> mesh.nodes[elem[i]], Val(4))
u_elem = ntuple(i -> Vec{3}((u[3*elem[i]-2], u[3*elem[i]-1], u[3*elem[i]])), Val(4))
# Element residual
r_elem = [zero(Vec{3,Float64}) for _ in 1:4]
assemble_element_residual!(r_elem, X, u_elem, material, element_states[elem_idx])
# Scatter to global
for (i, node) in enumerate(elem)
r[3*node-2] += r_elem[i][1]
r[3*node-1] += r_elem[i][2]
r[3*node] += r_elem[i][3]
end
end
return r
end
# ============================================================================
# Matrix-Free Operator for GMRES
# ============================================================================
"""
Matrix-free Jacobian-vector product: J·v (R(u+ε·v) - R(u))/ε
"""
struct MatrixFreeJacobian
u::Vector{Float64} # Current solution
r::Vector{Float64} # Current residual R(u)
mesh::Any # Mesh data
material::Any # Material
element_states::Vector # Plastic states (one vector per element)
ε::Float64 # Finite difference step
# Temporary storage
u_pert::Vector{Float64}
r_pert::Vector{Float64}
end
function MatrixFreeJacobian(u, r, mesh, material, element_states, ε=1e-7)
u_pert = similar(u)
r_pert = similar(r)
return MatrixFreeJacobian(u, r, mesh, material, element_states, ε, u_pert, r_pert)
end
"""
Apply J·v using finite differences.
"""
function apply_jacobian!(result, J::MatrixFreeJacobian, v)
# Perturbed solution: u + ε·v
@. J.u_pert = J.u + J.ε * v
# Need to copy states for perturbation (don't modify original)
states_pert = [copy(states) for states in J.element_states]
# Perturbed residual: R(u + ε·v)
assemble_global_residual!(J.r_pert, J.u_pert, J.mesh, J.material, states_pert)
# Finite difference: (R(u+ε·v) - R(u))/ε
@. result = (J.r_pert - J.r) / J.ε
return result
end
# ============================================================================
# GMRES Solver (Matrix-Free)
# ============================================================================
"""
Simple GMRES implementation (matrix-free).
Solves J·du = -r for Newton correction.
"""
function gmres_solve!(du, J::MatrixFreeJacobian, r, max_iter=50, tol=1e-6)
n = length(r)
# Krylov subspace basis
V = [zeros(n) for _ in 1:max_iter+1]
H = zeros(max_iter + 1, max_iter)
# Initial residual
fill!(du, 0.0)
V[1] .= -r # We're solving J·du = -r
β = norm(V[1])
V[1] ./= β
# Givens rotations storage
g = zeros(max_iter + 1)
g[1] = β
c = zeros(max_iter)
s = zeros(max_iter)
for j in 1:max_iter
# Arnoldi: w = J·v_j
w = zeros(n)
apply_jacobian!(w, J, V[j])
# Modified Gram-Schmidt
for i in 1:j
H[i, j] = dot(w, V[i])
w .-= H[i, j] .* V[i]
end
H[j+1, j] = norm(w)
if H[j+1, j] > 1e-14
V[j+1] .= w ./ H[j+1, j]
end
# Apply previous Givens rotations
for i in 1:j-1
temp = c[i] * H[i, j] + s[i] * H[i+1, j]
H[i+1, j] = -s[i] * H[i, j] + c[i] * H[i+1, j]
H[i, j] = temp
end
# Compute new Givens rotation
ρ = sqrt(H[j, j]^2 + H[j+1, j]^2)
c[j] = H[j, j] / ρ
s[j] = H[j+1, j] / ρ
H[j, j] = ρ
H[j+1, j] = 0.0
# Update residual norm
g[j+1] = -s[j] * g[j]
g[j] = c[j] * g[j]
residual_norm = abs(g[j+1])
if residual_norm < tol * β
# Back-solve upper triangular system
y = zeros(j)
for i in j:-1:1
y[i] = g[i]
for k in i+1:j
y[i] -= H[i, k] * y[k]
end
y[i] /= H[i, i]
end
# Form solution: du = V * y
for i in 1:j
du .+= y[i] .* V[i]
end
@printf(" GMRES converged in %d iterations (res: %.2e)\n", j, residual_norm)
return true
end
end
@printf(" GMRES did NOT converge after %d iterations\n", max_iter)
return false
end
# ============================================================================
# Anderson Acceleration
# ============================================================================
"""
Anderson acceleration for Newton iterations.
Mixes previous iterates to accelerate convergence.
"""
mutable struct AndersonAccelerator
m::Int # Mixing depth
X::Vector{Vector{Float64}} # Previous iterates
F::Vector{Vector{Float64}} # Previous residuals
iter::Int # Current iteration
end
function AndersonAccelerator(n::Int, m::Int=5)
X = [zeros(n) for _ in 1:m]
F = [zeros(n) for _ in 1:m]
return AndersonAccelerator(m, X, F, 0)
end
"""
Apply Anderson mixing to compute next iterate.
"""
function anderson_step!(acc::AndersonAccelerator, u_new, f_new)
acc.iter += 1
if acc.iter == 1
# First iteration - no mixing
return copy(u_new)
end
# Number of previous iterates to use
k = min(acc.iter - 1, acc.m)
# Store current iterate
idx = mod1(acc.iter, acc.m)
acc.X[idx] .= u_new
acc.F[idx] .= f_new
if k == 1
# Not enough history - just return current
return copy(u_new)
end
# Build ΔF matrix (differences of residuals)
ΔF = zeros(length(f_new), k - 1)
for i in 1:k-1
idx_curr = mod1(acc.iter - i + 1, acc.m)
idx_prev = mod1(acc.iter - i, acc.m)
ΔF[:, i] .= acc.F[idx_curr] .- acc.F[idx_prev]
end
# Solve least-squares: min ||ΔF·θ - f_new||
θ = ΔF \ f_new
# Mixed iterate: u = u_new - Σ θ_i (u_{k-i} - u_{k-i-1})
u_mixed = copy(u_new)
for i in 1:k-1
idx_curr = mod1(acc.iter - i + 1, acc.m)
idx_prev = mod1(acc.iter - i, acc.m)
u_mixed .-= θ[i] .* (acc.X[idx_curr] .- acc.X[idx_prev])
end
return u_mixed
end
# ============================================================================
# Newton Solver with Anderson Acceleration
# ============================================================================
"""
Newton solver with GMRES and Anderson acceleration.
This is the COMPLETE PIPELINE.
"""
function solve_newton_krylov_anderson!(u, mesh, material, element_states;
max_iter=20, tol=1e-6, anderson_depth=5)
println("\n" * "="^70)
println("Newton-Krylov-Anderson Solver")
println("="^70)
n = length(u)
r = zeros(n)
du = zeros(n)
# Anderson accelerator
anderson = AndersonAccelerator(n, anderson_depth)
for iter in 1:max_iter
# Assemble residual
assemble_global_residual!(r, u, mesh, material, element_states)
r_norm = norm(r)
@printf("Newton iter %2d: ||r|| = %.6e\n", iter, r_norm)
if r_norm < tol
println("✅ Converged!")
return true
end
# Matrix-free Jacobian
J = MatrixFreeJacobian(u, r, mesh, material, element_states)
# GMRES solve: J·du = -r
gmres_solve!(du, J, r)
# Line search parameter (simple version)
α = 1.0
u_new = u .+ α .* du
# Anderson acceleration (mix with previous iterates)
if anderson_depth > 0
u_new = anderson_step!(anderson, u_new, r)
end
# Update solution
u .= u_new
end
println("❌ Did NOT converge after $max_iter iterations")
return false
end
# ============================================================================
# Test Problem
# ============================================================================
# Simple mesh structure
struct SimpleMesh
nodes::Vector{Vec{3,Float64}}
elements::Vector{NTuple{4,Int}}
end
function create_test_mesh()
# Single Tet4 element
nodes = [
Vec{3}((0.0, 0.0, 0.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))
]
elements = [(1, 2, 3, 4)]
return SimpleMesh(nodes, elements)
end
function main()
println("\n" * "="^70)
println("Newton-Krylov-Anderson Reference Implementation (CPU)")
println("="^70)
# Mesh
mesh = create_test_mesh()
n_nodes = length(mesh.nodes)
n_dofs = 3 * n_nodes
# Material (perfect plasticity)
material = VonMisesPlasticity(
200e9, # E = 200 GPa
0.3, # ν = 0.3
250e6 # σ_y = 250 MPa
)
println("\n📦 Problem Setup:")
println(" Elements: $(length(mesh.elements))")
println(" Nodes: $n_nodes")
println(" DOFs: $n_dofs")
println(" Material: VonMises plasticity (σ_y = $(material.σ_y/1e6) MPa)")
# Initial guess (small displacement to trigger plasticity)
u = zeros(n_dofs)
u[4] = 0.002 # Node 2, x-direction (2mm - should exceed elastic limit)
# Initialize plastic states (one vector per element, one state per gauss point)
element_states = [[PlasticState() for _ in 1:length(GAUSS_TET4_1PT)]
for _ in 1:length(mesh.elements)]
# Solve
converged = solve_newton_krylov_anderson!(u, mesh, material, element_states,
max_iter=20, tol=1e-6, anderson_depth=3)
println("\n📊 Final Results:")
println(" Converged: $converged")
println(" u (first 6 DOFs): $(u[1:6])")
# Check plastic state
for (elem_idx, states) in enumerate(element_states)
for (gp_idx, state) in enumerate(states)
if state.α > 1e-10
@printf(" Element %d, GP %d: Plastic (α = %.6e)\n",
elem_idx, gp_idx, state.α)
end
end
end
println("\n✅ CPU REFERENCE COMPLETE!")
println("="^70)
end
# Run
main()
-387
View File
@@ -1,387 +0,0 @@
"""
Nodal Assembly CPU Reference Implementation
This implements the TWO-PHASE nodal assembly approach that will be ported to GPU:
Phase 1: Compute integration point data (stresses, material states)
Phase 2: Nodal assembly (matrix-free, no atomics)
Uses Tensors.jl throughout for natural tensor operations.
"""
using Tensors
using LinearAlgebra
using Printf
# Material state for plasticity
mutable struct PlasticState
ε_p::SymmetricTensor{2,3,Float64,6} # Plastic strain tensor
α::Float64 # Accumulated plastic strain
end
# Material properties
struct Material
E::Float64 # Young's modulus
ν::Float64 # Poisson's ratio
σ_y::Float64 # Yield stress
end
# Node-to-elements connectivity (CSR format)
struct NodeToElementsMap
ptr::Vector{Int} # Length: n_nodes + 1
data::Vector{Int} # Length: total connections
end
"""
Build CSR map: which elements touch each node?
"""
function build_node_to_elems(elements::Vector{NTuple{4,Int}}, n_nodes::Int)
# Count connections per node
counts = zeros(Int, n_nodes)
for elem in elements
for node in elem
counts[node] += 1
end
end
# Build CSR structure
ptr = cumsum([1; counts])
data = Vector{Int}(undef, sum(counts))
# Fill data array
offset = copy(ptr[1:end-1])
for (elem_idx, elem) in enumerate(elements)
for node in elem
data[offset[node]] = elem_idx
offset[node] += 1
end
end
return NodeToElementsMap(ptr, data)
end
"""
Tet4 shape function derivatives in reference coordinates (constant!)
"""
function tet4_shape_derivatives()
return (
Vec{3}((-1.0, -1.0, -1.0)), # dN1/dξ
Vec{3}((1.0, 0.0, 0.0)), # dN2/dξ
Vec{3}((0.0, 1.0, 0.0)), # dN3/dξ
Vec{3}((0.0, 0.0, 1.0)) # dN4/dξ
)
end
"""
Return mapping for von Mises perfect plasticity (using Tensors.jl!)
"""
function return_mapping_tensor(ε_total::SymmetricTensor{2,3,T},
state_old::PlasticState,
mat::Material) where T
# Elastic strain
ε_e = ε_total - state_old.ε_p
# Elastic predictor
λ = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
μ = mat.E / (2(1 + mat.ν))
I = one(ε_e)
σ_trial = λ * tr(ε_e) * I + 2μ * ε_e
# Deviatoric stress
σ_dev = dev(σ_trial)
σ_eq = sqrt(3 / 2 * (σ_dev σ_dev))
# Yield function
f = σ_eq - mat.σ_y
if f <= 0.0
# Elastic
return (σ_trial, state_old)
else
# Plastic - radial return
Δγ = f / (3μ)
n = σ_dev / σ_eq
σ = σ_trial - 2μ * Δγ * n
# Update plastic state
Δε_p = Δγ * n
ε_p_new = state_old.ε_p + Δε_p
α_new = state_old.α + Δγ
state_new = PlasticState(ε_p_new, α_new)
return (σ, state_new)
end
end
"""
Phase 1: Compute integration point data (stresses and material states)
One "thread" per integration point (in CPU version, just a loop)
"""
function compute_gp_data!(
σ_gp::Vector{SymmetricTensor{2,3,Float64,6}},
states_new::Vector{PlasticState},
u::Vector{Float64},
nodes::Matrix{Float64}, # Shape: 3 × n_nodes
elements::Vector{NTuple{4,Int}},
states_old::Vector{PlasticState},
mat::Material
)
n_elems = length(elements)
n_gps_per_elem = 4 # 4 Gauss points for Tet4
dN_dxi = tet4_shape_derivatives()
for elem_idx in 1:n_elems
# Extract element nodes
n1, n2, n3, n4 = elements[elem_idx]
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
u1 = Vec{3}((u[3*n1-2], u[3*n1-1], u[3*n1]))
u2 = Vec{3}((u[3*n2-2], u[3*n2-1], u[3*n2]))
u3 = Vec{3}((u[3*n3-2], u[3*n3-1], u[3*n3]))
u4 = Vec{3}((u[3*n4-2], u[3*n4-1], u[3*n4]))
# Jacobian (using tensor products!)
J = dN_dxi[1] X1 + dN_dxi[2] X2 + dN_dxi[3] X3 + dN_dxi[4] X4
invJ = inv(J)
# Physical derivatives
dN1_dx = invJ dN_dxi[1]
dN2_dx = invJ dN_dxi[2]
dN3_dx = invJ dN_dxi[3]
dN4_dx = invJ dN_dxi[4]
# Loop over Gauss points (for Tet4, same strain at all GPs since linear)
# In real code, would have different GP locations
for local_gp in 1:n_gps_per_elem
gp_idx = (elem_idx - 1) * n_gps_per_elem + local_gp
# Strain (using tensor products!)
ε = symmetric(dN1_dx u1 + dN2_dx u2 + dN3_dx u3 + dN4_dx u4)
# Material state update
state_old = states_old[gp_idx]
σ, state_new = return_mapping_tensor(ε, state_old, mat)
# Store results
σ_gp[gp_idx] = σ
states_new[gp_idx] = state_new
end
end
end
"""
Phase 2: Nodal assembly (matrix-free, no atomics!)
One "thread" per node (in CPU version, just a loop)
"""
function nodal_assembly!(
r::Vector{Float64},
σ_gp::Vector{SymmetricTensor{2,3,Float64,6}},
nodes::Matrix{Float64},
elements::Vector{NTuple{4,Int}},
node_to_elems::NodeToElementsMap
)
n_nodes = size(nodes, 2)
n_gps_per_elem = 4
dN_dxi = tet4_shape_derivatives()
# Gauss weights for Tet4 (standard 4-point quadrature)
gauss_weights = (1 / 24, 1 / 24, 1 / 24, 1 / 24)
fill!(r, 0.0)
for node_idx in 1:n_nodes
f_node = zero(Vec{3,Float64})
# Get elements touching this node (CSR traversal)
elem_start = node_to_elems.ptr[node_idx]
elem_end = node_to_elems.ptr[node_idx+1] - 1
# Loop over touching elements
for elem_offset in elem_start:elem_end
elem_idx = node_to_elems.data[elem_offset]
elem_nodes = elements[elem_idx]
# Find local node index in element
local_node = findfirst(==(node_idx), elem_nodes)
@assert local_node !== nothing "Node not found in element!"
# Recompute geometry (matrix-free approach!)
n1, n2, n3, n4 = elem_nodes
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
J = dN_dxi[1] X1 + dN_dxi[2] X2 + dN_dxi[3] X3 + dN_dxi[4] X4
detJ = det(J)
invJ = inv(J)
# Physical derivative for this node
dN_dx = invJ dN_dxi[local_node]
# Loop over Gauss points
for local_gp in 1:n_gps_per_elem
gp_idx = (elem_idx - 1) * n_gps_per_elem + local_gp
# Get stress at this GP
σ = σ_gp[gp_idx]
# Gauss weight
w = gauss_weights[local_gp]
# Accumulate force (using tensor contraction!)
f_node += (dN_dx σ) * (w * detJ)
end
end
# Write result (in GPU version, no atomics needed - this node is ours!)
r[3*node_idx-2] = f_node[1]
r[3*node_idx-1] = f_node[2]
r[3*node_idx] = f_node[3]
end
end
"""
Complete residual computation (two-phase approach)
"""
function compute_residual!(
r::Vector{Float64},
u::Vector{Float64},
nodes::Matrix{Float64},
elements::Vector{NTuple{4,Int}},
states_old::Vector{PlasticState},
mat::Material,
node_to_elems::NodeToElementsMap
)
n_gp = length(states_old)
# Storage for integration point data
σ_gp = Vector{SymmetricTensor{2,3,Float64,6}}(undef, n_gp)
states_new = Vector{PlasticState}(undef, n_gp)
# Phase 1: Compute integration point data
compute_gp_data!(σ_gp, states_new, u, nodes, elements, states_old, mat)
# Phase 2: Nodal assembly
nodal_assembly!(r, σ_gp, nodes, elements, node_to_elems)
return r, states_new
end
# ============================================================================
# Test Setup
# ============================================================================
function main()
println("\n" * "="^70)
println("Nodal Assembly CPU Reference - Using Tensors.jl")
println("="^70)
# Single Tet4 element
nodes = Float64[
0.0 1.0 0.0 0.0; # X coordinates
0.0 0.0 1.0 0.0; # Y coordinates
0.0 0.0 0.0 1.0 # Z coordinates
]
elements = [(1, 2, 3, 4)]
n_nodes = 4
n_elems = 1
n_gps = n_elems * 4 # 4 GPs per Tet4
# Material
mat = Material(
210e3, # E = 210 GPa (steel)
0.3, # ν = 0.3
250.0 # σ_y = 250 MPa
)
# Displacement (apply tension)
u = zeros(12)
u[4] = 0.01 # Move node 2 in X-direction
# Initial states (all elastic)
states_old = [PlasticState(zero(SymmetricTensor{2,3,Float64}), 0.0) for _ in 1:n_gps]
# Build node-to-elements map
println("\nBuilding node-to-elements map (CSR format)...")
node_to_elems = build_node_to_elems(elements, n_nodes)
println("CSR ptr: ", node_to_elems.ptr)
println("CSR data: ", node_to_elems.data)
# Verify each node touches exactly 1 element
for node_idx in 1:n_nodes
elem_start = node_to_elems.ptr[node_idx]
elem_end = node_to_elems.ptr[node_idx+1] - 1
n_touching = elem_end - elem_start + 1
touching_elems = node_to_elems.data[elem_start:elem_end]
println("Node $node_idx touches $n_touching element(s): $touching_elems")
end
# Compute residual (two-phase approach)
println("\n" * "-"^70)
println("Computing residual (two-phase nodal assembly)...")
println("-"^70)
r = zeros(12)
r, states_new = compute_residual!(r, u, nodes, elements, states_old, mat, node_to_elems)
println("\nResidual vector (internal forces):")
for i in 1:n_nodes
rx = r[3*i-2]
ry = r[3*i-1]
rz = r[3*i]
@printf("Node %d: [%12.6e, %12.6e, %12.6e]\n", i, rx, ry, rz)
end
println("\nResidual norm: ", norm(r))
# Check material states
println("\n" * "-"^70)
println("Material States at Gauss Points:")
println("-"^70)
for (gp_idx, state) in enumerate(states_new)
elem_idx = (gp_idx - 1) ÷ 4 + 1
local_gp = (gp_idx - 1) % 4 + 1
status = state.α > 0.0 ? "Plastic" : "Elastic"
@printf("Elem %d, GP %d: %s (α = %.6e)\n", elem_idx, local_gp, status, state.α)
end
# Test force balance (should sum to zero for internal forces)
println("\n" * "-"^70)
println("Force Balance Check:")
println("-"^70)
f_total = sum(reshape(r, 3, :), dims=2)
@printf("Sum of forces: [%.6e, %.6e, %.6e]\n", f_total[1], f_total[2], f_total[3])
@printf("Should be ≈ zero for internal forces (tol: 1e-10)\n")
if norm(f_total) < 1e-10
println("✅ Force balance: PASSED")
else
println("❌ Force balance: FAILED")
end
println("\n" * "="^70)
println("✅ Nodal assembly CPU reference complete!")
println("="^70)
println("\nNext step: Port this to GPU with CUDA.jl")
println(" Phase 1: @cuda compute_gp_data_kernel!(...)")
println(" Phase 2: @cuda nodal_assembly_kernel!(...)")
println("="^70 * "\n")
end
main()
-442
View File
@@ -1,442 +0,0 @@
"""
Nodal Assembly GPU Implementation
This is the GPU port of demos/nodal_assembly_cpu.jl using CUDA.jl.
TWO-PHASE APPROACH:
1. compute_gp_data_kernel!() - Compute integration point stresses and material states
2. nodal_assembly_kernel!() - Assemble residual at nodes (matrix-free, no atomics!)
Uses Tensors.jl throughout on GPU (CuArray{SymmetricTensor} works!)
"""
using CUDA
using Tensors
using LinearAlgebra
using Printf
# Material state for plasticity (GPU-compatible!)
struct PlasticState
ε_p::SymmetricTensor{2,3,Float64,6} # Plastic strain tensor
α::Float64 # Accumulated plastic strain
end
# Material properties
struct Material
E::Float64 # Young's modulus
ν::Float64 # Poisson's ratio
σ_y::Float64 # Yield stress
end
# Node-to-elements connectivity (CSR format)
struct NodeToElementsMap
ptr::CuArray{Int32,1} # Length: n_nodes + 1
data::CuArray{Int32,1} # Length: total connections
end
"""
Build CSR map: which elements touch each node?
"""
function build_node_to_elems_gpu(elements::Vector{NTuple{4,Int}}, n_nodes::Int)
# Count connections per node
counts = zeros(Int, n_nodes)
for elem in elements
for node in elem
counts[node] += 1
end
end
# Build CSR structure
ptr = cumsum([1; counts])
data = Vector{Int32}(undef, sum(counts))
# Fill data array
offset = copy(ptr[1:end-1])
for (elem_idx, elem) in enumerate(elements)
for node in elem
data[offset[node]] = elem_idx
offset[node] += 1
end
end
return NodeToElementsMap(CuArray(Int32.(ptr)), CuArray(data))
end
"""
Return mapping for von Mises perfect plasticity (using Tensors.jl on GPU!)
This function works identically on CPU and GPU!
"""
@inline function return_mapping_tensor(ε_total::SymmetricTensor{2,3,T},
state_old::PlasticState,
E, ν, σ_y) where T
# Elastic strain
ε_e = ε_total - state_old.ε_p
# Elastic predictor
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
I = one(ε_e)
σ_trial = λ * tr(ε_e) * I + 2μ * ε_e
# Deviatoric stress
σ_dev = dev(σ_trial)
σ_eq = sqrt(3 / 2 * (σ_dev σ_dev))
# Yield function
f = σ_eq - σ_y
if f <= T(0.0)
# Elastic
return (σ_trial, state_old)
else
# Plastic - radial return
Δγ = f / (3μ)
n = σ_dev / σ_eq
σ = σ_trial - 2μ * Δγ * n
# Update plastic state
Δε_p = Δγ * n
ε_p_new = state_old.ε_p + Δε_p
α_new = state_old.α + Δγ
state_new = PlasticState(ε_p_new, α_new)
return (σ, state_new)
end
end
"""
PHASE 1 GPU KERNEL: Compute integration point data
One thread per integration point!
"""
function compute_gp_data_kernel!(
σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1},
states_new::CuDeviceArray{PlasticState,1},
u::CuDeviceArray{Float64,1},
nodes::CuDeviceArray{Float64,2}, # Shape: 3 × n_nodes
elements::CuDeviceArray{Int32,2}, # Shape: 4 × n_elems
states_old::CuDeviceArray{PlasticState,1},
E, ν, σ_y
)
gp_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if gp_idx <= length(σ_gp)
# Map GP to element and local GP
elem_idx = (gp_idx - 1) ÷ 4 + 1 # 4 GPs per Tet4
# local_gp = (gp_idx - 1) % 4 + 1 # Not used yet (all GPs same for linear Tet4)
# Extract element nodes
n1 = elements[1, elem_idx]
n2 = elements[2, elem_idx]
n3 = elements[3, elem_idx]
n4 = elements[4, elem_idx]
# Node coordinates (using Tensors.jl Vec!)
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
# Displacements (using Tensors.jl Vec!)
u1 = Vec{3}((u[3*n1-2], u[3*n1-1], u[3*n1]))
u2 = Vec{3}((u[3*n2-2], u[3*n2-1], u[3*n2]))
u3 = Vec{3}((u[3*n3-2], u[3*n3-1], u[3*n3]))
u4 = Vec{3}((u[3*n4-2], u[3*n4-1], u[3*n4]))
# Shape derivatives (constant for Tet4)
dN1_dxi = Vec{3}((-1.0, -1.0, -1.0))
dN2_dxi = Vec{3}((1.0, 0.0, 0.0))
dN3_dxi = Vec{3}((0.0, 1.0, 0.0))
dN4_dxi = Vec{3}((0.0, 0.0, 1.0))
# Jacobian (using tensor products!)
J = dN1_dxi X1 + dN2_dxi X2 + dN3_dxi X3 + dN4_dxi X4
invJ = inv(J)
# Physical derivatives (using tensor contractions!)
dN1_dx = invJ dN1_dxi
dN2_dx = invJ dN2_dxi
dN3_dx = invJ dN3_dxi
dN4_dx = invJ dN4_dxi
# Strain (using tensor products and symmetric!)
ε = symmetric(dN1_dx u1 + dN2_dx u2 + dN3_dx u3 + dN4_dx u4)
# Material state update (using Tensors.jl - works on GPU!)
state_old = states_old[gp_idx]
σ, state_new = return_mapping_tensor(ε, state_old, E, ν, σ_y)
# Store results
σ_gp[gp_idx] = σ
states_new[gp_idx] = state_new
end
return nothing
end
"""
PHASE 2 GPU KERNEL: Nodal assembly (matrix-free, no atomics!)
One thread per node!
"""
function nodal_assembly_kernel!(
r::CuDeviceArray{Float64,1},
σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1},
nodes::CuDeviceArray{Float64,2},
elements::CuDeviceArray{Int32,2},
node_to_elems_ptr::CuDeviceArray{Int32,1},
node_to_elems_data::CuDeviceArray{Int32,1}
)
node_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if node_idx <= size(nodes, 2)
# Accumulate forces from all elements touching this node
f_node = zero(Vec{3,Float64})
# Gauss weights for Tet4 (standard 4-point quadrature)
gauss_weight = 1.0 / 24.0
# Shape derivatives in reference coordinates
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))
)
# Get element range for this node (CSR format)
elem_start = node_to_elems_ptr[node_idx]
elem_end = node_to_elems_ptr[node_idx+1] - 1
# Loop over touching elements
for elem_offset in elem_start:elem_end
elem_idx = node_to_elems_data[elem_offset]
# Extract element nodes
n1 = elements[1, elem_idx]
n2 = elements[2, elem_idx]
n3 = elements[3, elem_idx]
n4 = elements[4, elem_idx]
# Find local node index in element
local_node = 1
if node_idx == n2
local_node = 2
elseif node_idx == n3
local_node = 3
elseif node_idx == n4
local_node = 4
end
# Recompute geometry (matrix-free!)
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
J = dN_dxi[1] X1 + dN_dxi[2] X2 + dN_dxi[3] X3 + dN_dxi[4] X4
detJ = det(J)
invJ = inv(J)
# Physical derivative for this node
dN_dx = invJ dN_dxi[local_node]
# Loop over Gauss points (4 per Tet4)
for local_gp in 1:4
gp_idx = (elem_idx - 1) * 4 + local_gp
# Get stress at this GP
σ = σ_gp[gp_idx]
# Accumulate force (using tensor contraction!)
f_node += (dN_dx σ) * (gauss_weight * detJ)
end
end
# Write result (no atomics - this node is ours!)
r[3*node_idx-2] = f_node[1]
r[3*node_idx-1] = f_node[2]
r[3*node_idx] = f_node[3]
end
return nothing
end
"""
Complete residual computation on GPU (two-phase approach)
"""
function compute_residual_gpu!(
r::CuArray{Float64,1},
u::CuArray{Float64,1},
nodes::CuArray{Float64,2},
elements::CuArray{Int32,2},
states_old::CuArray{PlasticState,1},
mat::Material,
node_to_elems::NodeToElementsMap
)
n_gp = length(states_old)
n_nodes = size(nodes, 2)
# Storage for integration point data (on GPU!)
σ_gp = CuArray{SymmetricTensor{2,3,Float64,6}}(undef, n_gp)
states_new = CuArray{PlasticState}(undef, n_gp)
# Phase 1: Compute integration point data
threads = 256
blocks = cld(n_gp, threads)
@cuda threads = threads blocks = blocks compute_gp_data_kernel!(
σ_gp, states_new,
u, nodes, elements,
states_old,
mat.E, mat.ν, mat.σ_y
)
# Phase 2: Nodal assembly
fill!(r, 0.0)
threads = 256
blocks = cld(n_nodes, threads)
@cuda threads = threads blocks = blocks nodal_assembly_kernel!(
r, σ_gp,
nodes, elements,
node_to_elems.ptr, node_to_elems.data
)
return r, states_new
end
# ============================================================================
# Test Setup
# ============================================================================
function main()
println("\n" * "="^70)
println("Nodal Assembly GPU Implementation - CUDA.jl + Tensors.jl")
println("="^70)
# Check CUDA availability
if !CUDA.functional()
println("❌ CUDA not available! This demo requires a GPU.")
return
end
println("\n✅ CUDA device: ", CUDA.name(CUDA.device()))
# Single Tet4 element
nodes_cpu = Float64[
0.0 1.0 0.0 0.0; # X coordinates
0.0 0.0 1.0 0.0; # Y coordinates
0.0 0.0 0.0 1.0 # Z coordinates
]
elements_cpu = [(1, 2, 3, 4)]
elements_mat = Int32[e[i] for i in 1:4, e in elements_cpu]
n_nodes = 4
n_elems = 1
n_gps = n_elems * 4 # 4 GPs per Tet4
# Material
mat = Material(
210e3, # E = 210 GPa (steel)
0.3, # ν = 0.3
250.0 # σ_y = 250 MPa
)
# Displacement (apply tension)
u_cpu = zeros(12)
u_cpu[4] = 0.01 # Move node 2 in X-direction
# Initial states (all elastic)
states_old_cpu = [PlasticState(zero(SymmetricTensor{2,3,Float64}), 0.0) for _ in 1:n_gps]
# Build node-to-elements map
println("\nBuilding node-to-elements map (CSR format)...")
node_to_elems = build_node_to_elems_gpu(elements_cpu, n_nodes)
# Transfer to GPU
println("Transferring data to GPU...")
nodes_gpu = CuArray(nodes_cpu)
elements_gpu = CuArray(elements_mat)
u_gpu = CuArray(u_cpu)
states_old_gpu = CuArray(states_old_cpu)
r_gpu = CUDA.zeros(Float64, 12)
# Compute residual on GPU
println("\n" * "-"^70)
println("Computing residual on GPU (two-phase nodal assembly)...")
println("-"^70)
r_gpu, states_new_gpu = compute_residual_gpu!(
r_gpu, u_gpu, nodes_gpu, elements_gpu,
states_old_gpu, mat, node_to_elems
)
# Transfer results back to CPU
r_cpu = Array(r_gpu)
states_new_cpu = Array(states_new_gpu)
println("\nResidual vector (internal forces):")
for i in 1:n_nodes
rx = r_cpu[3*i-2]
ry = r_cpu[3*i-1]
rz = r_cpu[3*i]
@printf("Node %d: [%12.6e, %12.6e, %12.6e]\n", i, rx, ry, rz)
end
println("\nResidual norm: ", norm(r_cpu))
# Check material states
println("\n" * "-"^70)
println("Material States at Gauss Points:")
println("-"^70)
for (gp_idx, state) in enumerate(states_new_cpu)
elem_idx = (gp_idx - 1) ÷ 4 + 1
local_gp = (gp_idx - 1) % 4 + 1
status = state.α > 0.0 ? "Plastic" : "Elastic"
@printf("Elem %d, GP %d: %s (α = %.6e)\n", elem_idx, local_gp, status, state.α)
end
# Test force balance
println("\n" * "-"^70)
println("Force Balance Check:")
println("-"^70)
f_total = sum(reshape(r_cpu, 3, :), dims=2)
@printf("Sum of forces: [%.6e, %.6e, %.6e]\n", f_total[1], f_total[2], f_total[3])
@printf("Should be ≈ zero for internal forces (tol: 1e-10)\n")
if norm(f_total) < 1e-10
println("✅ Force balance: PASSED")
else
println("❌ Force balance: FAILED")
end
# Compare with CPU reference
println("\n" * "="^70)
println("Comparing with CPU reference (demos/nodal_assembly_cpu.jl)...")
println("="^70)
println("\nExpected residual norms should match!")
println(" CPU reference: 727.208...")
println(" GPU result: ", norm(r_cpu))
println("\n" * "="^70)
println("✅ GPU nodal assembly complete!")
println("="^70)
println("\nNext steps:")
println(" 1. Benchmark GPU vs CPU performance")
println(" 2. Scale to realistic mesh sizes (10K+ elements)")
println(" 3. Integrate with Newton-Krylov solver")
println(" 4. Add line search for convergence")
println(" 5. Add preconditioning (Chebyshev-Jacobi → GMG)")
println("="^70 * "\n")
end
main()
-157
View File
@@ -1,157 +0,0 @@
"""
Test Tet10 Assembly on CPU First
=================================
Validate the JuliaFEM pattern (loop through shape functions, no B-matrix)
before moving to GPU.
"""
using LinearAlgebra
using Tensors
# Material
struct LinearElastic
E::Float64
ν::Float64
end
λ(mat::LinearElastic) = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2mat.ν))
μ(mat::LinearElastic) = mat.E / (2(1 + mat.ν))
function compute_stress_3d(material::LinearElastic, eps::SymmetricTensor{2,3,T}) where T
lambda_val = T(λ(material))
mu_val = T(μ(material))
I = one(eps)
sigma = lambda_val * tr(eps) * I + 2 * mu_val * eps
return sigma
end
# Tet10 Gauss quadrature
const GAUSS_TET4 = [
(Vec{3}((0.5854101966249685, 0.1381966011250105, 0.1381966011250105)), 0.25),
(Vec{3}((0.1381966011250105, 0.5854101966249685, 0.1381966011250105)), 0.25),
(Vec{3}((0.1381966011250105, 0.1381966011250105, 0.5854101966249685)), 0.25),
(Vec{3}((0.1381966011250105, 0.1381966011250105, 0.1381966011250105)), 0.25)
]
function tet10_shape_derivatives(xi, eta, zeta)::NTuple{10,Vec{3,Float64}}
lambda = 1 - xi - eta - zeta
# Vertex nodes
dN1 = Vec{3}((4 * lambda - 1, 4 * lambda - 1, 4 * lambda - 1))
dN2 = Vec{3}((4 * xi - 1, 0.0, 0.0))
dN3 = Vec{3}((0.0, 4 * eta - 1, 0.0))
dN4 = Vec{3}((0.0, 0.0, 4 * zeta - 1))
# Edge midpoints
dN5 = Vec{3}((4 * (1 - 2 * xi - eta - zeta), -4 * xi, -4 * xi))
dN6 = Vec{3}((4 * eta, 4 * xi, 0.0))
dN7 = Vec{3}((-4 * eta, 4 * (1 - xi - 2 * eta - zeta), -4 * eta))
dN8 = Vec{3}((-4 * zeta, -4 * zeta, 4 * (1 - xi - eta - 2 * zeta)))
dN9 = Vec{3}((4 * zeta, 0.0, 4 * xi))
dN10 = Vec{3}((0.0, 4 * zeta, 4 * eta))
return (dN1, dN2, dN3, dN4, dN5, dN6, dN7, dN8, dN9, dN10)
end
function compute_jacobian_tet10(dN_dxi, X)
# J = Σ dN_i ⊗ X_i (tensor products!)
return sum(dN_dxi[i] X[i] for i in 1:10)
end
# Compute strain using tensor products
function compute_strain_from_displacements(dN_dx, u)
# ∇u = Σ dN_i ⊗ u_i
gradu = sum(dN_dx[i] u[i] for i in 1:10)
# ε = sym(∇u)
return symmetric(gradu)
end
# Compute forces using tensors
function compute_nodal_forces_from_stress(dN_dx, sigma)
# f_i = dN_i · σ
return ntuple(i -> dN_dx[i] sigma, Val(10))
end
# Test
function main()
println("\n" * "="^70)
println("Tet10 Assembly Test (CPU)")
println("="^70)
# Single Tet10 element
X = (
Vec{3}((0.0, 0.0, 0.0)), # 1
Vec{3}((1.0, 0.0, 0.0)), # 2
Vec{3}((0.0, 1.0, 0.0)), # 3
Vec{3}((0.0, 0.0, 1.0)), # 4
Vec{3}((0.5, 0.0, 0.0)), # 5
Vec{3}((0.5, 0.5, 0.0)), # 6
Vec{3}((0.0, 0.5, 0.0)), # 7
Vec{3}((0.0, 0.0, 0.5)), # 8
Vec{3}((0.5, 0.0, 0.5)), # 9
Vec{3}((0.0, 0.5, 0.5)) # 10
)
# Displacements: Small perturbation
u = (
Vec{3}((0.0, 0.0, 0.0)), # Node 1
Vec{3}((0.001, 0.0, 0.0)), # Node 2 (1mm in x)
Vec{3}((0.0, 0.0, 0.0)), # Node 3
Vec{3}((0.0, 0.0, 0.0)), # Node 4
Vec{3}((0.0005, 0.0, 0.0)), # Node 5
Vec{3}((0.0005, 0.0, 0.0)), # Node 6
Vec{3}((0.0, 0.0, 0.0)), # Node 7
Vec{3}((0.0, 0.0, 0.0)), # Node 8
Vec{3}((0.0005, 0.0, 0.0)), # Node 9
Vec{3}((0.0, 0.0, 0.0)) # Node 10
)
material = LinearElastic(200e9, 0.3)
println("✅ Material: E=$(material.E/1e9) GPa, ν=$(material.ν)")
println("✅ Pattern: Tensors.jl tensor products (⊗, ⋅, sym)")
# Assemble element residual
r_elem = [zero(Vec{3}) for _ in 1:10]
for (xivec, w) in GAUSS_TET4
xi, eta, zeta = xivec[1], xivec[2], xivec[3]
# Shape function derivatives
dN_dxi = tet10_shape_derivatives(xi, eta, zeta)
# Jacobian
J = compute_jacobian_tet10(dN_dxi, X)
detJ = det(J)
invJ = inv(J)
# Physical derivatives: dN/dx = invJ · dN/dxi (tensor contraction)
dN_dx = ntuple(i -> invJ Vec{3}(dN_dxi[i]), Val(10))
# Compute strain: ε = sym(∇u) where ∇u = Σ dN_i ⊗ u_i
eps = compute_strain_from_displacements(dN_dx, u)
# Compute stress
sigma = compute_stress_3d(material, eps)
# Compute forces: f_i = dN_i · σ
f_contrib = compute_nodal_forces_from_stress(dN_dx, sigma)
# Accumulate
for i in 1:10
r_elem[i] += f_contrib[i] * (w * detJ)
end
end
r_total = vcat([r_elem[i][j] for i in 1:10 for j in 1:3]...)
println("\n📊 Results:")
println(" ||r||: $(norm(r_total))")
println(" r[1:6] (node 1-2, x,y,z): $(r_total[1:6])")
println("\n✅ CPU TEST COMPLETE!")
println("="^70 * "\n")
end
main()