diff --git a/demos/GPU_KERNEL_PLAN.md b/demos/GPU_KERNEL_PLAN.md deleted file mode 100644 index f48e2ff..0000000 --- a/demos/GPU_KERNEL_PLAN.md +++ /dev/null @@ -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!** diff --git a/demos/README_GPU_MPI.md b/demos/README_GPU_MPI.md deleted file mode 100644 index a73028d..0000000 --- a/demos/README_GPU_MPI.md +++ /dev/null @@ -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. diff --git a/demos/README_GPU_POC.md b/demos/README_GPU_POC.md deleted file mode 100644 index 6fe3995..0000000 --- a/demos/README_GPU_POC.md +++ /dev/null @@ -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. diff --git a/demos/README_KRYLOV_DEMO.md b/demos/README_KRYLOV_DEMO.md deleted file mode 100644 index 483961e..0000000 --- a/demos/README_KRYLOV_DEMO.md +++ /dev/null @@ -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. diff --git a/demos/README_TENSORS_CORRECTION.md b/demos/README_TENSORS_CORRECTION.md deleted file mode 100644 index 8aaf351..0000000 --- a/demos/README_TENSORS_CORRECTION.md +++ /dev/null @@ -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. diff --git a/demos/assembly_comparison_simple.jl b/demos/assembly_comparison_simple.jl deleted file mode 100644 index caf59a7..0000000 --- a/demos/assembly_comparison_simple.jl +++ /dev/null @@ -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) diff --git a/demos/cantilever_beam_demo.jl b/demos/cantilever_beam_demo.jl deleted file mode 100644 index 36f4c2d..0000000 --- a/demos/cantilever_beam_demo.jl +++ /dev/null @@ -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() diff --git a/demos/cantilever_beam_simple.jl b/demos/cantilever_beam_simple.jl deleted file mode 100644 index ee24255..0000000 --- a/demos/cantilever_beam_simple.jl +++ /dev/null @@ -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) diff --git a/demos/cantilever_cpu_comparison.jl b/demos/cantilever_cpu_comparison.jl deleted file mode 100644 index f5023fd..0000000 --- a/demos/cantilever_cpu_comparison.jl +++ /dev/null @@ -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) diff --git a/demos/cantilever_gmsh_gpu.jl b/demos/cantilever_gmsh_gpu.jl deleted file mode 100644 index a47e053..0000000 --- a/demos/cantilever_gmsh_gpu.jl +++ /dev/null @@ -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)") diff --git a/demos/cantilever_physics_gpu.jl b/demos/cantilever_physics_gpu.jl deleted file mode 100644 index 7cfdce5..0000000 --- a/demos/cantilever_physics_gpu.jl +++ /dev/null @@ -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") diff --git a/demos/gpu_assembly_poc.jl b/demos/gpu_assembly_poc.jl deleted file mode 100644 index e93ca86..0000000 --- a/demos/gpu_assembly_poc.jl +++ /dev/null @@ -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 diff --git a/demos/gpu_assembly_poc_tensors.jl b/demos/gpu_assembly_poc_tensors.jl deleted file mode 100644 index deb7f52..0000000 --- a/demos/gpu_assembly_poc_tensors.jl +++ /dev/null @@ -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 diff --git a/demos/gpu_assembly_tet10.jl b/demos/gpu_assembly_tet10.jl deleted file mode 100644 index ee0004f..0000000 --- a/demos/gpu_assembly_tet10.jl +++ /dev/null @@ -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 diff --git a/demos/gpu_elementset_demo.jl b/demos/gpu_elementset_demo.jl deleted file mode 100755 index d646a23..0000000 --- a/demos/gpu_elementset_demo.jl +++ /dev/null @@ -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)") diff --git a/demos/gpu_mpi_demo.jl b/demos/gpu_mpi_demo.jl deleted file mode 100644 index 73d654a..0000000 --- a/demos/gpu_mpi_demo.jl +++ /dev/null @@ -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() diff --git a/demos/gpu_mpi_mock.jl b/demos/gpu_mpi_mock.jl deleted file mode 100644 index 19f8a38..0000000 --- a/demos/gpu_mpi_mock.jl +++ /dev/null @@ -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() diff --git a/demos/gpu_nodal_assembly_demo.jl b/demos/gpu_nodal_assembly_demo.jl deleted file mode 100755 index 603f013..0000000 --- a/demos/gpu_nodal_assembly_demo.jl +++ /dev/null @@ -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() diff --git a/demos/gpu_only_demo.jl b/demos/gpu_only_demo.jl deleted file mode 100644 index 49962d3..0000000 --- a/demos/gpu_only_demo.jl +++ /dev/null @@ -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 diff --git a/demos/krylov_mpi_gpu_demo.jl b/demos/krylov_mpi_gpu_demo.jl deleted file mode 100644 index 2bdddbd..0000000 --- a/demos/krylov_mpi_gpu_demo.jl +++ /dev/null @@ -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() diff --git a/demos/newton_krylov_anderson_cpu.jl b/demos/newton_krylov_anderson_cpu.jl deleted file mode 100644 index 70cb483..0000000 --- a/demos/newton_krylov_anderson_cpu.jl +++ /dev/null @@ -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() diff --git a/demos/nodal_assembly_cpu.jl b/demos/nodal_assembly_cpu.jl deleted file mode 100644 index 8e95f89..0000000 --- a/demos/nodal_assembly_cpu.jl +++ /dev/null @@ -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() diff --git a/demos/nodal_assembly_gpu.jl b/demos/nodal_assembly_gpu.jl deleted file mode 100644 index f0e5b57..0000000 --- a/demos/nodal_assembly_gpu.jl +++ /dev/null @@ -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() diff --git a/demos/test_tet10_cpu.jl b/demos/test_tet10_cpu.jl deleted file mode 100644 index f1d62f3..0000000 --- a/demos/test_tet10_cpu.jl +++ /dev/null @@ -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()