diff --git a/docs/blog/immutability_tldr.md b/docs/blog/immutability_tldr.md deleted file mode 100644 index be7e4b2..0000000 --- a/docs/blog/immutability_tldr.md +++ /dev/null @@ -1,139 +0,0 @@ ---- -title: "Copying is Faster Than Mutating: A Counterintuitive Performance Win" -author: "Jukka Aho" -date: "2025-11-09" -categories: ["Performance", "Benchmarks"] -tags: ["immutability", "type-stability", "quick-reference"] -description: "TL;DR version: How immutable elements are 130x faster with zero allocations" ---- - -## TL;DR - -We made our FEM code **130x faster** by making it immutable. Yes, copying everything is faster than mutating in place. No, we're not crazy. We have benchmarks. - -## The "Stupid" Idea - -**Old code (mutable Dict):** - -```julia -element.fields[:E] = 210e9 # Mutate in place - "fast" -``` - -**New code (immutable NamedTuple):** - -```julia -element = update(element, E=210e9) # Copy entire struct - "slow" -``` - -Which do you think is faster? - -## The Shocking Results - -```text -Field access: 41x faster (immutable) -Assembly loop: 130x faster (immutable) -1000 element mesh: 120x faster (immutable) - -Memory allocations: 70,000 → 0 (immutable) -``` - -**Immutable is 100x faster AND uses zero memory.** - -## The Secret: Type Stability - -```julia -# Dict{Symbol,Any} - Type unstable -element.fields[:E] # Compiler: "What type is this? 🤷" -# Cost: hash lookup + pointer chase + runtime dispatch ≈ 45 nanoseconds - -# NamedTuple{(:E,:ν),Tuple{Float64,Float64}} - Type stable -element.fields.E # Compiler: "Float64 at offset 0. Got it." -# Cost: inline to CPU register ≈ 1 nanosecond -``` - -**45x slower just to read a field.** Multiply by millions of accesses in FEM assembly. - -## The Compiler Magic - -When you write: - -```julia -element = ImmutableElement((1,2,3,4), (E=210e9, ν=0.3)) -element = update(element, temperature=293.15) -``` - -The compiler sees: - -- Old element not used → reuse stack space -- New element same size → copy is one assignment -- All types known → inline everything -- Result: **Zero heap allocations, SIMD vectorization, GPU-ready** - -When you write: - -```julia -element.fields[:temperature] = 293.15 -``` - -The Dict must: - -- Compute hash of `:temperature` -- Check if key exists (pointer chasing) -- Maybe resize Dict (heap allocation) -- Store as `Any` → runtime dispatch on next access -- Result: **Heap allocations, type instability, CPU-only** - -## Real-World Impact - -**Assemble 10,000 element mesh:** - -| Implementation | Time | Memory | GPU | -|----------------|------|---------|-----| -| Dict (mutable) | 2.4s | 450 MB, 7M allocs | ✗ | -| NamedTuple (immutable) | **0.02s** | **0 MB, 0 allocs** | ✓ | - -Interactive vs coffee break. Million-element mesh vs out-of-memory. GPU vs CPU-only. - -## The Lesson - -Your programming intuition is from 1990s C/C++: - -- ✓ Mutation is fast ← **TRUE IN C** -- ✓ Copying is slow ← **TRUE IN C** -- ✗ Type doesn't matter ← **FALSE IN MODERN COMPILERS** - -2025 reality: - -- **Type stability is everything** -- Compiler optimizes away struct copies -- Mutation breaks type inference -- Immutability enables GPU acceleration - -## Try It Yourself - -```bash -git clone https://github.com/JuliaFEM/JuliaFEM.jl -cd JuliaFEM.jl -julia benchmarks/element_immutability_benchmark.jl -``` - -Full article: `docs/blog/immutability_performance.md` - -## Bottom Line - -We made the "wrong" choice (copy everything, mutate nothing) and got: - -- 130x faster code -- Zero allocations -- GPU compatibility -- Better parallelization - -**Copying > Mutating. Immutability > Mutation. Type stability > Everything.** - -Measure, don't assume. The evidence is in the benchmarks. - ---- - -*JuliaFEM 1.0 architecture, November 2025* -*Benchmark: Intel i7-12700K, Julia 1.12.1* -*Full results in repository* diff --git a/docs/blog/krylov_nodal_assembly.jl b/docs/blog/krylov_nodal_assembly.jl deleted file mode 100644 index e497e75..0000000 --- a/docs/blog/krylov_nodal_assembly.jl +++ /dev/null @@ -1,415 +0,0 @@ -# # Krylov Subspace Iterations Meet Nodal Assembly: A Revolution in Contact Mechanics -# -# **Author:** Jukka Aho -# **Date:** November 2025 -# **Status:** Vision and demonstration of JuliaFEM v1.0 architecture - -# ## The Clever Combination Nobody Talks About -# -# Here's something that should be obvious but isn't: **Krylov subspace methods combined -# with nodal assembly are the natural way to solve contact problems**. Yet almost every -# FEM code does it the hard way—assembling global matrices element-by-element and then -# complaining about memory usage. -# -# Why is nodal assembly + Krylov so brilliant for contact mechanics? -# -# 1. **Contact is inherently nodal**: When you write down the weak form of contact, -# the constraints appear at nodes, not elements. Contact forces, gaps, friction—all -# defined node-to-node. -# -# 2. **Krylov methods don't need the matrix**: They only need the matrix-vector product -# operator. You never have to form the global matrix. Just give me `y = A*x` and I'm -# happy. -# -# 3. **Nodal assembly gives you the rows**: Building matrix rows node-by-node is the -# most natural way to incorporate nodal contact constraints. No scatter/gather -# gymnastics needed. -# -# The result? **O(N) memory instead of O(N²), and contact constraints that fall out -# naturally from the formulation.** -# -# Traditional FEM codes can't do this because they're locked into element assembly -# paradigms from the 1970s. We're not. - -# ## A Broader Vision: Nodal Material Modeling -# -# Here's where it gets controversial. Today, everyone computes material state -# (stress, plastic strain, damage) at **integration points**. This feels natural -# because that's where we evaluate integrals, right? -# -# **Wrong. It's backwards.** -# -# Think about it: Why should the material state depend on the numerical integration -# scheme? You can't get analytical solutions at the element level because you've -# *assumed* you'll use Gaussian quadrature. The physics is now **constrained by the -# numerical method**. That's insane! -# -# ### The Nodal Material State Hypothesis -# -# **I claim that material modeling should also be nodal**, for the same reasons contact -# is nodal: -# -# 1. **Variational consistency**: The weak form naturally places material response at -# nodes when you do it properly. Integration points are an implementation detail. -# -# 2. **Physical meaning**: Nodes represent physical points in space. Integration points? -# They're mathematical constructs that change when you pick a different quadrature rule. -# -# 3. **Scalability**: Nodal material state scales linearly with problem size. -# Integration point state scales with elements × points per element. -# -# 4. **Contact-material coupling**: When contact happens, material state at the contact -# node matters. Why store it somewhere else and interpolate? -# -# ### Why Nobody Believes This (Yet) -# -# Every material scientist will tell you I'm crazy. "You need integration points for -# plasticity!" "What about locking?" "This violates the patch test!" -# -# **I will show them they're wrong.** Not today, but it's coming. The math works out -# when you do the variational formulation correctly. It just requires thinking beyond -# 1970s element technology. -# -# For now, we focus on contact (where nodal is already accepted), and we build the -# infrastructure that will eventually support nodal materials too. - -# ## The Krylov Advantage: Solving Unsymmetric Systems -# -# Here's another key insight: **Real problems are unsymmetric**. -# -# - **Material nonlinearity**: Tangent stiffness from plasticity is usually unsymmetric -# - **Contact**: Contact contributions are inherently unsymmetric (one-sided constraints) -# - **Large deformations**: Geometric nonlinearity introduces unsymmetry -# -# Traditional FEM codes use direct solvers (LU decomposition) which don't care about -# symmetry but scale as O(N³). Iterative solvers designed for symmetric problems -# (Conjugate Gradient) fail on unsymmetric systems. -# -# **Enter GMRES**: Generalized Minimal Residual method. It solves unsymmetric systems -# as long as they're invertible (positive definite is enough). Combined with nodal -# assembly, you get: -# -# - O(N·iter) time complexity (vs O(N³) for direct) -# - O(N) memory (vs O(N²) for storing full matrix) -# - Handles unsymmetry naturally -# - Works with contact, plasticity, large deformation—everything - -# ## Demonstration: GMRES on Unsymmetric System -# -# Let's prove this works with a simple example: 10×10 positive definite but -# unsymmetric system, solved with GMRES using nodal assembly pattern. - -using LinearAlgebra -using Random -using Printf - -println("="^70) -println("GMRES + Nodal Assembly: Unsymmetric System Demo") -println("="^70) -println() - -# ### Problem Setup -# -# Create a positive definite but unsymmetric matrix. This mimics what you get from -# contact mechanics or material nonlinearity. - -Random.seed!(42) -N = 10 - -# Start with symmetric positive definite -A_sym = rand(N, N) -A_sym = A_sym' * A_sym + 10.0 * I(N) - -# Add small unsymmetric part (mimics contact or material nonlinearity) -# Keep it small to maintain positive definiteness -A_unsym = rand(N, N) * 0.1 -A = A_sym + A_unsym - -# Check if positive definite (all eigenvalues positive and real) -evals = eigvals(A) -evals_real = real.(evals) -all_real = all(abs.(imag.(evals)) .< 1e-10) -all_positive = all(evals_real .> 0) - -println("Matrix properties:") -println(" Size: $(N)×$(N)") -println(" Symmetric: ", issymmetric(A)) -if all_real - println(" Eigenvalues (real): ", evals_real) - println(" All positive: ", all_positive) -else - println(" Eigenvalues: ", evals) - println(" All positive: ", all_positive) -end -println(" Condition number: ", cond(A)) -println() - -# ### Exact Solution -# -# We know the answer—this lets us verify convergence. - -x_exact = Float64[i for i in 1:N] -b = A * x_exact - -println("Exact solution: x = [1, 2, 3, ..., $N]") -println() - -# ### Nodal Assembly Pattern -# -# Define the row-by-row assembly interface. In real FEM, `get_row(i)` would -# assemble contributions from all elements connected to node `i`. - -""" - get_row(A, i) -> Vector{Float64} - -Nodal assembly: return the i-th row of the system matrix. -In real FEM, this would sum contributions from all elements touching node i. -""" -function get_row(A::Matrix{Float64}, i::Int) - return A[i, :] -end - -""" - matvec_nodal(A, x) -> Vector{Float64} - -Matrix-vector product using nodal assembly. -Computes y = A*x by assembling and using one row at a time. -""" -function matvec_nodal(A::Matrix{Float64}, x::Vector{Float64}) - n = length(x) - y = zeros(Float64, n) - - for i in 1:n - row = get_row(A, i) - y[i] = dot(row, x) - end - - return y -end - -# Test the nodal matvec -x_test = ones(N) -y_test = matvec_nodal(A, x_test) -y_direct = A * x_test -println("Nodal matvec test:") -println(" Error vs direct: ", norm(y_test - y_direct)) -println(" ✓ Nodal assembly working correctly") -println() - -# ### GMRES Implementation -# -# Simplified GMRES for demonstration. Production code would use Krylov.jl, -# but this shows the core algorithm clearly. - -""" - gmres_simple(A, b, x0; maxiter=100, tol=1e-10) - -Simplified GMRES using nodal assembly pattern. -Only needs matrix-vector product—never forms full matrix. -""" -function gmres_simple(A::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64}; - maxiter::Int=100, tol::Float64=1e-10) - n = length(b) - x = copy(x0) - - # Arnoldi iteration vectors - V = zeros(Float64, n, maxiter + 1) - H = zeros(Float64, maxiter + 1, maxiter) - - # Initial residual - r = b - matvec_nodal(A, x) - β = norm(r) - V[:, 1] = r / β - - # Store residual history - residuals = Float64[β] - - println("GMRES iteration:") - @printf(" Initial residual: %.6e\n", β) - println() - - for j in 1:maxiter - # Arnoldi: build orthonormal basis for Krylov subspace - w = matvec_nodal(A, V[:, j]) - - # Modified Gram-Schmidt orthogonalization - 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 - - # Solve least squares problem: min ||β*e₁ - H*y|| - e1 = zeros(j + 1) - e1[1] = β - - # Use QR factorization (simple, stable) - Hj = H[1:j+1, 1:j] - y = Hj \ e1 - - # Update solution - x_new = x0 + V[:, 1:j] * y - - # Compute residual - r = b - matvec_nodal(A, x_new) - res_norm = norm(r) - push!(residuals, res_norm) - - reduction = 100.0 * (1.0 - res_norm / β) - @printf(" Iteration %3d: residual = %.6e (reduction: %.2f%%)\n", - j, res_norm, reduction) - - if res_norm < tol - println() - println(" ✓ Converged in $j iterations") - return x_new, j, res_norm, residuals - end - - x = x_new - end - - println() - println(" ⚠ Did not converge in $maxiter iterations") - return x, maxiter, norm(b - matvec_nodal(A, x)), residuals -end - -# ### Solve with GMRES - -x0 = zeros(N) -x_solution, iters, final_res, res_history = gmres_simple(A, b, x0, maxiter=50, tol=1e-10) - -println() - -# ### Verification - -error_abs = norm(x_solution - x_exact) -error_rel = error_abs / norm(x_exact) - -println("="^70) -println("Verification Results") -println("="^70) -println() -println("Solution comparison:") -println(" Exact: ", join([@sprintf("%.3f", x) for x in x_exact], ", ")) -println(" Computed: ", join([@sprintf("%.3f", x) for x in x_solution], ", ")) -println() -println("Error metrics:") -@printf(" Absolute error: %.6e\n", error_abs) -@printf(" Relative error: %.6e\n", error_rel) -@printf(" Final residual: %.6e\n", final_res) -println(" Iterations: $iters") -println() - -if error_rel < 1e-6 - println("✅ VERIFICATION PASSED") -else - println("❌ VERIFICATION FAILED") -end -println() - -# ## Key Insights from This Demonstration - -println("="^70) -println("Why This Matters for JuliaFEM v1.0") -println("="^70) -println() - -println(""" -1. **Unsymmetric systems are solved naturally** - - Matrix is positive definite but unsymmetric ✓ - - GMRES converges in $iters iterations ✓ - - Solution accurate to 1e-$(Int(round(-log10(error_rel)))) relative error ✓ - - This is what real contact/plasticity problems look like - -2. **Nodal assembly works perfectly** - - Never formed global matrix explicitly - - Only used get_row(i) interface—one row at a time - - Memory: O(N) instead of O(N²) - - In real FEM: get_row(i) assembles from elements touching node i - -3. **Krylov methods scale** - - This demo: $N×$N system, $iters iterations - - Scales to millions: 1M×1M system, ~100 iterations typical - - Time: O(N·iter) vs O(N³) for direct solvers - - Memory: O(N) vs O(N²) for storing full matrix - -4. **Contact mechanics fits naturally** - - Contact constraints modify rows for contact nodes - - No special treatment needed—just part of get_row(i) - - Nodal formulation, nodal assembly, nodal constraints - - Everything at the same level—beautiful! - -5. **Foundation for future: Nodal materials** - - Same infrastructure supports nodal material state - - Material history at nodes, not integration points - - Physically meaningful, numerically efficient - - Controversial today, obvious tomorrow -""") - -println("="^70) -println() - -# ## The Path Forward -# -# This demonstration proves the concept works. For JuliaFEM v1.0: -# -# ### Immediate (Months 1-3) -# - Implement `get_row(node_id, elements)` for real element assembly -# - Integrate Krylov.jl for production-quality GMRES -# - Add preconditioning (Jacobi, ILU) for faster convergence -# - Handle contact constraints in row modification -# -# ### Near-term (Months 4-6) -# - Matrix-free operators with GPU acceleration -# - Distributed assembly across MPI ranks -# - Strong scaling studies (speedup vs number of processes) -# - Contact mechanics validation (Hertz, patch tests) -# -# ### Long-term (Months 7-12) -# - Nodal material state experiments (plasticity at nodes) -# - Compare integration-point vs nodal material models -# - Publish results showing nodal materials work -# - Prove the material scientists wrong 😎 -# -# ### Vision (Beyond v1.0) -# - Complete nodal formulation: geometry, contact, materials -# - Demonstrate 10M DOF contact problems on multi-GPU clusters -# - Show that element-centric thinking was 20th century -# - Lead the field into 21st century FEM - -# ## Conclusion -# -# **Krylov subspace iterations + nodal assembly is not just a technical choice—it's -# the philosophically correct approach to contact mechanics.** -# -# Contact is nodal. Constraints are nodal. Solution method should be nodal. -# -# Traditional FEM uses element assembly because that's how it was done in 1970 -# (before iterative solvers were practical). We're not constrained by history. -# -# **And soon we'll show that materials should be nodal too.** The math works. -# The numerics work (as shown in this demo). The physics makes sense. -# -# It just requires thinking clearly about what the weak form actually says, -# rather than cargo-culting element assembly from outdated textbooks. -# -# *Welcome to JuliaFEM v1.0. Where we assemble by nodes, solve with Krylov, -# and refuse to be constrained by integration point theology.* -# -# --- -# -# **References:** -# - Saad, Y. (2003). *Iterative Methods for Sparse Linear Systems*. SIAM. -# - Wriggers, P. (2006). *Computational Contact Mechanics*. Springer. -# - Aho, J. (2025). "Why Material Scientists Are Wrong About Integration Points" -# (forthcoming, controversy expected 😉) - -println("Demo complete. For production code, see:") -println(" - demos/krylov_mpi_gpu_demo.jl (distributed multi-GPU solver)") -println(" - docs/book/nodal_assembly_multigpu.md (strategic document)") -println() diff --git a/docs/design/IMMUTABILITY.md b/docs/design/IMMUTABILITY.md deleted file mode 100644 index dcd3de9..0000000 --- a/docs/design/IMMUTABILITY.md +++ /dev/null @@ -1,491 +0,0 @@ ---- -title: "Element Immutability: Design Decision and Rationale" -author: "Jukka Aho" -date: "2025-11-09" -status: "IMPLEMENTED" -phase: "Phase 1B" -categories: ["Design", "Architecture"] -tags: ["immutability", "performance", "design-decision"] -benchmark: "benchmarks/element_immutability_benchmark.jl" ---- - -## Executive Summary - -JuliaFEM 1.0 adopts **immutable elements with type-stable fields** as a core architectural decision. While this appears counterintuitive (requiring element copies instead of in-place mutation), benchmarks demonstrate **40-130x performance improvement** over the mutable Dict-based approach. - -**Key Results:** - -- Field access: **40x faster** (1ns vs 45ns per read) -- Assembly loop: **130x faster** (9ns vs 1,124ns per element) -- Large mesh: **120x faster** (0.01ms vs 1.2ms for 1000 elements) -- Memory: **Zero allocations** in hot path (vs 70,000 allocations) -- GPU/HPC: **Compatible** (all bits types vs pointers) - ---- - -## The Counterintuitive API Change - -### Old API (Mutable, Dict-based) - -```julia -# Create element with mutable fields -element = Element(Tet10, [1,2,3,4,5,6,7,8,9,10]) - -# Add fields dynamically -update!(element, "E", 210e9) -update!(element, "ν", 0.3) -update!(element, "temperature", 293.15) - -# Fields stored in Dict{Symbol,Any} - type unstable! -element.fields # → Dict(:E => 210e9, :ν => 0.3, :temperature => 293.15) -``` - -**Pros:** Familiar, flexible, feels efficient (no copies) -**Cons:** Type-unstable, 100ns Dict lookup overhead, no GPU support - -### New API (Immutable, Type-stable) - -```julia -# Create element with type-stable fields -element = Element(Lagrange{Tetrahedron,2}, (1,2,3,4,5,6,7,8,9,10), - fields=(E=210e9, ν=0.3)) - -# Update returns NEW element (immutable) -element = update(element, temperature=293.15) - -# Fields stored in NamedTuple - type stable! -element.fields # → (E=210e9, ν=0.3, temperature=293.15) -typeof(element.fields) # → NamedTuple{(:E,:ν,:temperature), Tuple{Float64,Float64,Float64}} -``` - -**Pros:** Type-stable, 1ns access, GPU-compatible, zero allocations -**Cons:** Requires element copy (but compiler optimizes away!) - ---- - -## Why Immutability Wins - -### 1. Type Stability is Everything - -In FEM assembly, field access happens **millions of times**: - -```julia -# Assembly loop: 10 integration points × 1000 elements = 10,000 field accesses -for element in mesh - for ip in integration_points - E = element.fields[:E] # Dict lookup: 45ns EACH TIME - ν = element.fields[:ν] # Another 45ns - # ... compute stiffness - end -end -``` - -**Mutable (Dict):** `45ns × 20,000 = 900µs` (Dict lookups) -**Immutable (Tuple):** `1ns × 20,000 = 20µs` (direct access) - -**Result:** 45x speedup just from field access! - -### 2. Compiler Optimizations - -Type-stable code enables: - -- **Inlining:** Field access becomes single instruction -- **SIMD:** Vectorization across multiple elements -- **Constant propagation:** Compiler knows exact types -- **Stack allocation:** No heap allocations for small structs - -Example: Assembly loop with immutable elements **completely inlines**: - -```julia -# Before optimization (conceptual): -E = element.fields.E # Field access -λ = E * ν / ... # Material computation - -# After optimization (actual machine code): -λ = 210e9 * 0.3 / ... # Constants folded, direct computation! -``` - -### 3. Zero Allocations - -**Mutable elements:** Every field update allocates - -```julia -julia> @benchmark update!(element, "temperature", 293.15) -Allocs: 100 # One allocation per update! -Memory: 1600 bytes -``` - -**Immutable elements:** Stack allocation only - -```julia -julia> @benchmark element = update(element, temperature=293.15) -Allocs: 0 # Compiler optimizes to stack! -Memory: 0 bytes -``` - -**Why?** Modern Julia compiler recognizes stack-only pattern and eliminates heap allocations entirely. - -### 4. GPU/HPC Compatibility - -**Mutable elements with Dict:** - -```julia -struct MutableElement - fields::Dict{Symbol,Any} # POINTER → cannot transfer to GPU -end -``` - -**Immutable elements with NamedTuple:** - -```julia -struct ImmutableElement{F} - fields::F # All bits types → can transfer to GPU! -end -``` - -GPU kernels require: - -- No pointers (CPU memory → GPU memory not allowed) -- No dynamic dispatch (GPU can't call CPU functions) -- All data as bits types (can be copied to GPU) - -Only immutable, type-stable elements satisfy these requirements. - ---- - -## Benchmark Results - -Run: `julia --project=. benchmarks/element_immutability_benchmark.jl` - -### Field Access (1000 reads) - -| Implementation | Time/read | Speedup | -|---------------|-----------|---------| -| Mutable (Dict) | 45ns | 1x (baseline) | -| Immutable (Tuple) | 1ns | **40x** | - -### Field Update (100 writes) - -| Implementation | Time/update | Allocations | -|---------------|-------------|-------------| -| Mutable (mutate) | 12ns | 100 | -| Immutable (copy) | 0.03ns | 0 | - -**Surprise:** Creating new structs is **400x faster** than mutating Dict! - -### Assembly Loop (single element) - -| Implementation | Time | Allocations | -|---------------|------|-------------| -| Mutable | 1,124ns | 69 | -| Immutable | 9ns | 0 | - -**Speedup:** **130x faster** - -### Large Mesh (1000 elements) - -| Implementation | Time | Memory | -|----------------|--------|--------| -| Mutable | 1.2ms | 1.1 MB | -| Immutable | 0.01ms | 0 KB | - -**Speedup:** **120x faster**, zero allocations - ---- - -## Common Misconceptions - -### "Copying structs is expensive" - -**False.** Small structs (< 128 bytes) are stack-allocated: - -```julia -# This looks like it copies: -new_element = update(old_element, temperature=300.0) - -# But actually compiles to: -# mov rax, [old_fields] # Load old fields -# mov [new_fields], rax # Store to new location (STACK!) -# mov [new_fields+24], 300.0 # Update temperature field -``` - -No heap allocation, no GC pressure, just register/stack operations. - -### "I need mutable fields for time integration" - -**False.** Time-varying fields should be stored separately: - -```julia -# Bad: Time history in element (mutable) -element.fields[:temperature] = [293.15, 300.0, 310.0] # Vector → allocates - -# Good: Time history separate (immutable element) -struct TimeHistory - times::Vector{Float64} - temperatures::Vector{Float64} -end - -element = Element(..., fields=(E=210e9, ν=0.3)) # Constant -history = TimeHistory([0.0, 1.0, 2.0], [293.15, 300.0, 310.0]) # Mutable separately -``` - -Element stays immutable (fast), history is mutable (when needed). - -### "Functional programming is slow" - -**False in Julia.** Persistent data structures (like Clojure) are slow because they allocate on heap. Julia's immutable structs are stack-allocated and get optimized away by compiler. - -```julia -# This code: -e1 = Element(..., fields=(E=210e9,)) -e2 = update(e1, ν=0.3) -e3 = update(e2, ρ=7850.0) - -# Compiles to: -# Stack allocation: -# [E] [ν] [ρ] -# 210e9 0.3 7850.0 ← Single struct on stack! -``` - ---- - -## Design Patterns - -### Pattern 1: Initialization with Fields - -```julia -# Create element with all known fields upfront -element = Element(Lagrange{Triangle,1}, (1,2,3), - fields=(E=210e9, ν=0.3, thickness=0.01)) -``` - -### Pattern 2: Progressive Updates - -```julia -# Start with minimal fields -element = Element(Lagrange{Triangle,1}, (1,2,3), fields=(E=210e9,)) - -# Add fields as computed (returns new element) -element = update(element, ν=0.3) -element = update(element, temperature=compute_temperature(element)) -``` - -### Pattern 3: Batch Updates - -```julia -# Update multiple fields at once (efficient!) -element = update(element, - temperature=300.0, - stress=(σ_xx=100e6, σ_yy=50e6, σ_xy=0.0), - plastic_strain=0.001) -``` - -### Pattern 4: Field Inheritance - -```julia -# Reuse fields from another element -base_fields = (E=210e9, ν=0.3, ρ=7850.0) - -elem1 = Element(Lagrange{Triangle,1}, (1,2,3), fields=base_fields) -elem2 = Element(Lagrange{Triangle,1}, (4,5,6), fields=base_fields) -# Both share same type → compiler can optimize across elements! -``` - -### Pattern 5: Conditional Fields - -```julia -# Different elements can have different field sets -function create_element(topology, conn, use_plasticity) - if use_plasticity - fields = (E=210e9, ν=0.3, yield_stress=250e6) - else - fields = (E=210e9, ν=0.3) - end - return Element(topology, conn, fields=fields) -end -``` - ---- - -## Migration Guide (Old → New) - -### Old Code (Mutable) - -```julia -# Create element -element = Element(Tet10, [1,2,3,4,5,6,7,8,9,10]) - -# Add fields -update!(element, "E", 210e9) -update!(element, "ν", 0.3) - -# Access fields -E = element.fields[:E] -``` - -### New Code (Immutable) - -```julia -# Create element with fields -element = Element(Lagrange{Tetrahedron,2}, (1,2,3,4,5,6,7,8,9,10), - fields=(E=210e9, ν=0.3)) - -# Update returns new element -element = update(element, temperature=293.15) - -# Access fields (type-stable!) -E = element.fields.E -``` - -### Key Changes - -1. **Creation:** Include fields at construction time -2. **Update:** Assign result: `element = update(element, ...)` -3. **Access:** Use dot syntax: `element.fields.E` not `element.fields[:E]` -4. **Types:** Prefer NamedTuple over Dict: `(E=210e9,)` not `Dict(:E => 210e9)` - ---- - -## Implementation Details - -### Element Definition - -```julia -struct Element{N,NIP,F,B} <: AbstractElement{F,B} - id::UInt - connectivity::NTuple{N,UInt} # Immutable tuple - integration_points::NTuple{NIP,IP} # Immutable tuple - fields::F # Type-stable! (NamedTuple or struct) - basis::B # Type-stable! -end -``` - -### Update Implementation - -```julia -function update(element::Element, new_fields::NamedTuple) - # Merge old and new fields - updated_fields = merge(element.fields, new_fields) - - # Create new element (same connectivity, new fields) - return Element{N,NIP,typeof(updated_fields),B}( - element.id, - element.connectivity, - element.integration_points, - updated_fields, - element.basis - ) -end - -# Convenience syntax -update(element; kwargs...) = update(element, values(kwargs)) -``` - -### Memory Layout - -```julia -# Old mutable element (heap): -MutableElement -├── id: UInt64 (8 bytes on stack) -├── connectivity: Vector (24 bytes pointer → heap) -└── fields: Dict (24 bytes pointer → heap) - ↓ - [Heap allocations] - -# New immutable element (stack): -ImmutableElement -├── id: UInt64 (8 bytes) -├── connectivity: Tuple (40 bytes, inline) -└── fields: NamedTuple (24 bytes, inline) - ├── E: Float64 (8 bytes) - ├── ν: Float64 (8 bytes) - └── ρ: Float64 (8 bytes) - -Total: 72 bytes, all on stack, cache-friendly! -``` - ---- - -## Future Work - -### Phase 2: Time-Varying Fields - -Currently, fields are static. For time integration: - -```julia -# Option 1: External time history (current approach) -struct TimeVaryingField{T} - times::Vector{Float64} - values::Vector{T} -end - -# Element stays immutable -element = Element(..., fields=(E=210e9,)) -temperature_history = TimeVaryingField([0.0, 1.0], [293.15, 300.0]) - -# Option 2: Functional fields (future) -element = Element(..., fields=( - E=210e9, - temperature=t -> 293.15 + 10.0*t # Function of time -)) -``` - -### Phase 3: GPU Kernels - -With immutable elements, GPU assembly becomes possible: - -```julia -using CUDA - -# Transfer elements to GPU (all bits types!) -d_elements = CuArray(elements) -d_nodes = CuArray(nodes) - -# GPU kernel (parallel over elements) -@cuda threads=256 blocks=ceil(Int, n_elements/256) assemble_kernel!( - d_K, d_elements, d_nodes -) - -# No CPU synchronization needed - immutable = no race conditions! -``` - -### Phase 4: SIMD Vectorization - -Type-stable elements enable SIMD: - -```julia -# Process 4 elements simultaneously (AVX2) -function assemble_batch(elements::NTuple{4,Element}) - @simd for i in 1:4 - E = elements[i].fields.E # Vectorized load! - # ... assembly computation - end -end -``` - ---- - -## Conclusion - -**Immutability is not a compromise - it's an optimization.** - -Key takeaways: - -1. **Type stability dominates performance** in tight loops -2. **Compiler optimizations** make immutability free -3. **Zero allocations** eliminate GC pressure -4. **GPU/HPC compatibility** requires immutability -5. **Functional patterns** are fast in Julia - -The 40-130x speedup speaks for itself. Immutable elements are the foundation for high-performance, GPU-ready FEM in JuliaFEM 1.0. - ---- - -## References - -- Benchmark: `benchmarks/element_immutability_benchmark.jl` -- Implementation: `src/elements/elements.jl` -- Discussion: GitHub Issue #XXX (TBD) -- Related: `docs/design/FIELDS_DESIGN.md` (Phase 3) - -**Last Updated:** November 9, 2025