From b9534e2c55f9755a0e194cb03cd714ea1f8b72e8 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 15 Dec 2025 16:19:59 +0200 Subject: [PATCH] test(solvers): add modal analysis new API visionary test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New 560-line test file for modal analysis NEW API (test-driven development): - Visionary tests for ModalSolver and EigenvalueProblem - Tests single Tet4 element modal analysis - Tests with and without geometric stiffness - Tests heat transfer modal analysis - Extensive documentation of intended API design - Tests currently skipped (@test_skip) until implementation complete - Documents eigenvalue problem K φ = λ M φ - Documents mass matrix assembly and mode shape normalization This test file serves as both test suite and API design documentation for the new modal analysis solver interface. --- test/solvers/test_modal_analysis_new_api.jl | 560 ++++++++++++++++++++ 1 file changed, 560 insertions(+) create mode 100644 test/solvers/test_modal_analysis_new_api.jl diff --git a/test/solvers/test_modal_analysis_new_api.jl b/test/solvers/test_modal_analysis_new_api.jl new file mode 100644 index 0000000..d7ba957 --- /dev/null +++ b/test/solvers/test_modal_analysis_new_api.jl @@ -0,0 +1,560 @@ +""" +# Modal Analysis - NEW API (Test-Driven Development) + +**What:** Shows how modal (eigenvalue) analysis SHOULD work with the NEW API + +**Why:** +- **Eigenvalue problems** - Different solver type (not Newton-Krylov) +- **Natural frequencies** - Validates dynamic behavior +- **Geometric stiffness** - Prestress effects on frequencies +- **Multiple physics** - Works for elasticity AND heat transfer +- **Mass matrix required** - First time we need M, not just K + +**NEW API Concepts:** +1. **ModalSolver** - Eigenvalue solver (not LinearSolver or NewtonSolver) +2. **EigenvalueProblem** - K φ = λ M φ (generalized eigenvalue problem) +3. **Geometric stiffness** - K_g from initial stress/displacement +4. **Modal shapes** - φ_i (eigenvectors = mode shapes) +5. **Natural frequencies** - ω_i = √λ_i + +**Test Problems:** + +## Test 1: Single Tet4 Element +- 4-node tetrahedron with 3 nodes fixed +- Validates basic eigenvalue computation +- Tests with and without geometric stiffness + +## Test 2: Heat Transfer Modal Analysis +- Two Quad4 elements (plane heat) +- Eigenvalues of thermal diffusion operator +- Validates modal analysis works for temperature field + +## Test 3: Vibrating Beam (Future) +- Fixed-fixed beam +- Compare natural frequencies to analytical solution +- ω_i = λ_i² √(EI/ρA) + +**Expected Behavior (when implemented):** +✅ Eigenvalues computed correctly +✅ Mode shapes orthogonal (φ_i^T M φ_j = δ_ij) +✅ Geometric stiffness affects frequencies +✅ Works for both elasticity and heat transfer +✅ Sparse eigenvalue solver (only need first k modes) + +**Status:** 🚧 VISIONARY TEST - Implementation in progress +""" + +using Test +using JuliaFEM +using Tensors +using LinearAlgebra +using Statistics + +@testset "Modal Analysis - NEW API (TDD)" begin + + # ============================================================================= + # SINGLE ELEMENT: TET4 WITHOUT GEOMETRIC STIFFNESS + # ============================================================================= + + @testset "Single Tet4 Modal Analysis (Visionary)" begin + @test_skip begin # Skip until implemented + + # Geometry + X = Dict( + 1 => Vec(2.0, 3.0, 4.0), + 2 => Vec(6.0, 3.0, 2.0), + 3 => Vec(2.0, 5.0, 1.0), + 4 => Vec(4.0, 3.0, 6.0) + ) + + # Initial displacement (for geometric stiffness later) + u0 = Dict( + 1 => Vec(0.0, 0.0, 0.0), + 2 => Vec(0.0, 0.0, 0.0), + 3 => Vec(0.0, 0.0, 0.0), + 4 => Vec(0.25, 0.25, 0.25) + ) + + # Material + material = LinearElastic(E=96.0, ν=1 / 3, ρ=420.0) + + # Physics + elastic_physics = ContinuumPhysics{Displacement}( + material=material, + formulation=FullThreeD(), + finite_strain=false, + geometric_stiffness=false # No K_g initially + ) + + # Create element + tet_element = create_element( + Tet4, (1, 2, 3, 4), + geometry=X, + displacement=u0 + ) + + # Domain + domain = Domain( + name="TET", + elements=[tet_element], + physics=elastic_physics + ) + + # Boundary conditions (3 nodes fixed) + bc_fixed = DirichletBC( + name="FIXED_FACE", + nodes=[1, 2, 3], + dof=:displacement, + values=[0.0, 0.0, 0.0] + ) + + # NEW: EigenvalueProblem (not static or transient) + problem = EigenvalueProblem( + domains=[domain], + boundary_conditions=[bc_fixed] + ) + + # NEW: ModalSolver + solver = ModalSolver( + n_modes=2, # Number of modes to compute + which=:SM, # Smallest Magnitude (or :LM for largest) + method=:arpack # ARPACK, KrylovKit, or direct + ) + + # Solve: K φ = λ M φ + solution = solve!(problem, solver) + + # Extract eigenvalues and eigenvectors + λ = solution.eigenvalues # [λ_1, λ_2] + φ = solution.eigenvectors # [φ_1, φ_2] + + # Validate eigenvalues (without geometric stiffness) + @test isapprox(λ, [4 / 3, 1 / 3], rtol=1e-6) + + # Validate orthogonality: φ_i^T M φ_j = δ_ij + for i in 1:2 + for j in 1:2 + orthogonality = φ[i]' * solution.mass_matrix * φ[j] + expected = (i == j) ? 1.0 : 0.0 + @test isapprox(orthogonality, expected, atol=1e-6) + end + end + + end + end + + # ============================================================================= + # GEOMETRIC STIFFNESS EFFECT + # ============================================================================= + + @testset "Tet4 with Geometric Stiffness (Visionary)" begin + @test_skip begin + + # Same setup as before, but enable geometric stiffness + elastic_physics = ContinuumPhysics{Displacement}( + material=LinearElastic(E=96.0, ν=1 / 3, ρ=420.0), + formulation=FullThreeD(), + finite_strain=false, + geometric_stiffness=true # NEW: Enable K_g + ) + + problem = EigenvalueProblem( + domains=[domain], + boundary_conditions=[bc_fixed], + initial_displacement=u0 # Required for K_g + ) + + solver = ModalSolver(n_modes=2, which=:SM) + solution = solve!(problem, solver) + + λ = solution.eigenvalues + + # With geometric stiffness: K_eff = K + K_g(u0) + # Eigenvalues should be different! + @test isapprox(λ, [5 / 3, 2 / 3], rtol=1e-6) + + # Geometric stiffness changes natural frequencies + @test λ[1] > 4 / 3 # Stiffened by prestress + @test λ[2] > 1 / 3 + + end + end + + # ============================================================================= + # HEAT TRANSFER MODAL ANALYSIS + # ============================================================================= + + @testset "Heat Transfer Modal (Visionary)" begin + @test_skip begin + + # Geometry: Two quad elements stacked vertically + X = Dict( + 1 => Vec(0.0, 0.0), + 2 => Vec(1.0, 0.0), + 3 => Vec(1.0, 3.0), + 4 => Vec(0.0, 3.0), + 5 => Vec(0.0, 3.0), + 6 => Vec(1.0, 3.0), + 7 => Vec(1.0, 9.0), + 8 => Vec(0.0, 9.0) + ) + + # Material (thermal) + thermal_material = ThermalMaterial( + conductivity=36.0, + density=6.0, + specific_heat=1.0 + ) + + # Physics + thermal_physics = ContinuumPhysics{Temperature}( + material=thermal_material, + formulation=PlaneHeat() + ) + + # Elements + elem1 = create_element(Quad4, (1, 2, 3, 4), geometry=X) + elem2 = create_element(Quad4, (4, 3, 7, 8), geometry=X) + + # Domain + domain = Domain( + name="THERMAL", + elements=[elem1, elem2], + physics=thermal_physics + ) + + # Boundary conditions (fixed temperature at ends) + bc_bottom = DirichletBC( + nodes=[1, 2], + dof=:temperature, + value=0.0 + ) + + bc_top = DirichletBC( + nodes=[7, 8], + dof=:temperature, + value=0.0 + ) + + # Eigenvalue problem for heat equation + # M ∂T/∂t + K T = 0 + # Modal: K φ = λ M φ + problem = EigenvalueProblem( + domains=[domain], + boundary_conditions=[bc_bottom, bc_top] + ) + + solver = ModalSolver(n_modes=1, which=:SM) + solution = solve!(problem, solver) + + λ = solution.eigenvalues + + # First eigenvalue should be 1.0 (analytical) + @test isapprox(λ[1], 1.0, rtol=1e-6) + + end + end + + # ============================================================================= + # NATURAL FREQUENCIES FROM EIGENVALUES + # ============================================================================= + + @testset "Natural Frequency Computation (Visionary)" begin + @test_skip begin + + # Solve modal problem + solution = solve_modal_problem() + + # Extract eigenvalues + λ = solution.eigenvalues # [rad²/s²] + + # Natural frequencies: ω = √λ [rad/s] + ω = sqrt.(λ) + + # Convert to Hz: f = ω / (2π) + f = ω ./ (2π) + + # Validate units and values + @test all(λ .≥ 0) # Eigenvalues non-negative + @test all(ω .≥ 0) # Frequencies non-negative + @test all(f .≥ 0) + + # For simple beam: f_1 ≈ 1-10 Hz (typical) + @test 1.0 < f[1] < 100.0 + + end + end + + # ============================================================================= + # MODE SHAPE VISUALIZATION (FUTURE) + # ============================================================================= + + @testset "Mode Shape Properties (Visionary)" begin + @test_skip begin + + solution = solve_modal_problem() + + φ = solution.eigenvectors # Mode shapes + M = solution.mass_matrix + K = solution.stiffness_matrix + λ = solution.eigenvalues + + # Property 1: Orthogonality w.r.t. mass matrix + # φ_i^T M φ_j = δ_ij + for i in 1:length(φ) + for j in 1:length(φ) + orth_M = φ[i]' * M * φ[j] + expected = (i == j) ? 1.0 : 0.0 + @test isapprox(orth_M, expected, atol=1e-6) + end + end + + # Property 2: Eigenvalue equation + # K φ_i = λ_i M φ_i + for i in 1:length(φ) + Kφ = K * φ[i] + λMφ = λ[i] * (M * φ[i]) + @test Kφ ≈ λMφ rtol = 1e-6 + end + + # Property 3: Mode shapes normalized + # φ_i^T M φ_i = 1 + for i in 1:length(φ) + norm_M = φ[i]' * M * φ[i] + @test isapprox(norm_M, 1.0, rtol=1e-6) + end + + end + end + + # ============================================================================= + # MASS MATRIX ASSEMBLY + # ============================================================================= + + @testset "Mass Matrix Assembly (Visionary)" begin + # Pseudo-code showing mass matrix assembly + + println("\n" * "="^70) + println("MASS MATRIX ASSEMBLY (NODAL APPROACH)") + println("="^70) + + mass_assembly_pseudo = """ + # Similar to stiffness, but uses density ρ and N^T N + for node_i in nodes + for elem in node_to_elements[node_i] + for node_j in elem.nodes + # Mass matrix block (3×3 for displacement) + M_ij = ∫_Ω ρ N_i N_j dΩ + + # For displacement (3D): + # M_ij = m_ij * I_3×3 (often lumped) + + # Consistent mass (full integration) + M_block = compute_mass_block(elem, node_i, node_j, ρ) + + # Or lumped mass (diagonal only) + if lumped + M_block = (i == j) ? diag(M_block) : 0 + end + + M_nodal[node_i, node_j] += M_block + end + end + end + """ + + println(mass_assembly_pseudo) + println("="^70) + println("✓ Mass matrix M = ∫_Ω ρ N^T N dΩ") + println("✓ Consistent mass: Full integration (more accurate)") + println("✓ Lumped mass: Diagonal only (faster, explicit dynamics)") + println("✓ Nodal assembly works same as stiffness!") + println("="^70) + end + + # ============================================================================= + # KEY ARCHITECTURAL INSIGHTS + # ============================================================================= + + println("\n" * "="^70) + println("MODAL ANALYSIS ARCHITECTURE INSIGHTS (NEW API)") + println("="^70) + println("✓ Modal analysis is EigenvalueProblem (K φ = λ M φ)") + println("✓ Requires both stiffness K AND mass M matrices") + println("✓ ModalSolver uses sparse eigenvalue methods (ARPACK)") + println("✓ Geometric stiffness K_g affects natural frequencies") + println("✓ Works for ANY physics (displacement, temperature, etc.)") + println("✓ Mode shapes φ orthogonal w.r.t. mass matrix") + println("✓ Natural frequencies ω = √λ, in Hz: f = ω/(2π)") + println("✓ Same nodal assembly pattern for M as for K!") + println("="^70) + +end + +""" +# IMPLEMENTATION NOTES + +## Generalized Eigenvalue Problem + +**Mathematical formulation:** + +K φ = λ M φ + +where: +- K = stiffness matrix (from elastic/thermal energy) +- M = mass matrix (from kinetic/thermal energy) +- φ = eigenvector (mode shape) +- λ = eigenvalue (ω² for vibrations, α for heat) + +**For structural vibrations:** +- ω = √λ (natural frequency in rad/s) +- f = ω/(2π) (natural frequency in Hz) +- φ = displacement mode shape + +**For heat diffusion:** +- α = λ (thermal diffusivity eigenvalue) +- φ = temperature mode shape + +## Stiffness Matrix + +**Elasticity:** + +K = ∫_Ω B^T C B dΩ + +where: +- B = strain-displacement matrix +- C = elasticity tensor + +**Heat transfer:** + +K = ∫_Ω k ∇N^T ∇N dΩ + +where: +- k = thermal conductivity + +## Mass Matrix + +**Elasticity (consistent):** + +M = ∫_Ω ρ N^T N dΩ + +where ρ = density. + +**Elasticity (lumped):** + +M_ii = ∑_elements ∫_Ω_e ρ N_i dΩ_e + +(Diagonal only, faster for explicit dynamics) + +**Heat transfer:** + +M = ∫_Ω ρc N^T N dΩ + +where: +- ρ = density +- c = specific heat + +## Geometric Stiffness + +**Definition:** Stiffness contribution from initial stress state. + +K_g = ∫_Ω G^T σ₀ G dΩ + +where: +- G = geometric matrix (relates δε to ∇(δu)) +- σ₀ = initial stress + +**Effect:** Prestress changes natural frequencies: +- Tension → higher frequencies (stiffening) +- Compression → lower frequencies (softening) + +**Total effective stiffness:** + +K_eff = K + K_g + +Then solve: K_eff φ = λ M φ + +## Sparse Eigenvalue Solvers + +**Problem:** For large systems, computing ALL eigenvalues is expensive. + +**Solution:** Sparse eigenvalue methods (compute only k << n eigenvalues). + +**Methods:** +1. **ARPACK** - Arnoldi iteration (standard in Julia) +2. **KrylovKit.jl** - Modern Krylov methods +3. **Lanczos** - For symmetric problems (K and M symmetric) + +**JuliaFEM approach:** + +```julia +using Arpack + +# Solve K φ = λ M φ for k smallest eigenvalues +λ, φ = eigs(K, M, nev=k, which=:SM) +``` + +## Mode Shape Normalization + +**Goal:** Normalize eigenvectors for convenience. + +**Mass normalization (standard):** + +φ_i^T M φ_i = 1 + +**Advantages:** +- Orthogonality: φ_i^T M φ_j = δ_ij +- Energy interpretation clear +- Modal damping ratios well-defined + +**Implementation:** + +```julia +function normalize_modes!(φ, M) + for i in 1:length(φ) + # Compute φ_i^T M φ_i + norm_M = φ[i]' * M * φ[i] + + # Normalize + φ[i] ./= sqrt(norm_M) + end +end +``` + +## Nodal Assembly for Mass Matrix + +```julia +function assemble_mass_nodal!(M_nodal, ρ, elements, node_to_elements) + Threads.@threads for node_i in 1:n_nodes + for elem_idx in node_to_elements[node_i] + elem = elements[elem_idx] + + for node_j in elem.nodes + # Mass matrix block (3×3 for displacement) + M_ij = compute_mass_block(elem, node_i, node_j, ρ) + + # Add to global (no atomic on diagonal if i = j) + if node_i == node_j + M_nodal[node_i, node_i] += M_ij # Direct write + else + atomic_add!(M_nodal[node_i, node_j], M_ij) + end + end + end + end +end +``` + +**Key:** SAME PATTERN as stiffness assembly! + +## Next Steps + +1. Implement `EigenvalueProblem` type +2. Implement `ModalSolver` with ARPACK integration +3. Implement mass matrix assembly (nodal) +4. Implement geometric stiffness computation +5. Add mode shape normalization +6. Validate against analytical solutions +7. Performance benchmarks (large systems) + +"""