From 05397e02007d9e864c284a1644213de23f460661 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:25:42 +0300 Subject: [PATCH] chore(test): remove crazy THM multifield experiment Delete exploratory thermo-hydro-mechanical DOF tests outside CI coverage. - Drop `test/dofs/test_thm_crazy_multifield.jl`. --- test/dofs/test_thm_crazy_multifield.jl | 399 ------------------------- 1 file changed, 399 deletions(-) delete mode 100644 test/dofs/test_thm_crazy_multifield.jl diff --git a/test/dofs/test_thm_crazy_multifield.jl b/test/dofs/test_thm_crazy_multifield.jl deleted file mode 100644 index d8c1148..0000000 --- a/test/dofs/test_thm_crazy_multifield.jl +++ /dev/null @@ -1,399 +0,0 @@ -""" -πŸš€πŸš€πŸš€ THE CRAZIEST DEMO EVER: Multi-Field THM-E Coupling πŸš€πŸš€πŸš€ - -THIS IS IT! The ULTIMATE demonstration of coupled multi-physics using the -NEW multi-field Element API! - -ONE element with FOUR field types on FOUR different entity types: -- Temperature (Float64) at VERTICES -- Displacement (Vec{3}) at VERTICES -- Pore pressure (Float64) at CELLS -- Electric potential (Float64) at EDGES - -Physics: Fully coupled THM-E system: -1. Heat equation: βˆ‚T/βˆ‚t - ΞΊΞ”T = Q (conduction) -2. Darcy flow: βˆ‡Β·q = 0, q = -k(βˆ‡p + ρg) (pore pressure) -3. Elasticity: -βˆ‡Β·Οƒ = f, Οƒ = C:(Ξ΅ - Ξ±_T*T - Ξ±_p*p) (thermal + pore expansion) -4. Electrokinetics: -βˆ‡Β·(Οƒ_eβˆ‡Ο†) = 0 (electric potential on edges) - -Use case: Geomechanics, soil consolidation, electrokinetic remediation, - nuclear waste storage, CO2 sequestration, geothermal energy - -πŸŽ‰ NEW MULTI-FIELD API: ------------------------ -ONE element creation call! -Natural field access: elem.dof_indices.T, .u, .p, .Ο† -Type-safe field names! -No manual DOF management! - -THIS IS THE FUTURE! πŸš€ -""" - -using JuliaFEM -using Test -using Tensors -using LinearAlgebra -using SparseArrays -using Printf - -@testset "πŸš€πŸš€πŸš€ CRAZIEST DEMO EVER: Multi-Field THM-E" begin - println("\n" * "="^70) - println("πŸš€πŸš€πŸš€ THE CRAZIEST DEMO EVER: MULTI-FIELD THM-E πŸš€πŸš€πŸš€") - println("="^70) - - # Create 3D mesh: Two tetrahedra forming a simple domain - # Tet 1: nodes (1, 2, 3, 4) - # Tet 2: nodes (2, 3, 4, 5) - shares face with tet 1 - nodes = [ - Vec{3,Float64}((0.0, 0.0, 0.0)), # Node 1 - Vec{3,Float64}((1.0, 0.0, 0.0)), # Node 2 - Vec{3,Float64}((0.5, 1.0, 0.0)), # Node 3 - Vec{3,Float64}((0.5, 0.5, 1.0)), # Node 4 - Vec{3,Float64}((1.5, 0.5, 0.5)), # Node 5 - ] - connectivity = [ - (UInt32(1), UInt32(2), UInt32(3), UInt32(4)), # Tet 1 - (UInt32(2), UInt32(3), UInt32(4), UInt32(5)), # Tet 2 - ] - mesh = Mesh{Tetrahedron{4}}(nodes, connectivity) - - println("\n3D Mesh: 2 tetrahedra") - println(" 5 nodes, 9 edges, 7 faces, 2 cells") - - println("\n" * "="^70) - println("πŸŽ‰ CREATING MULTI-FIELD ELEMENTS - THE NEW WAY!") - println("="^70) - - # Define multi-field specification - field_spec = NamedTuple{(:T, :u, :p, :Ο†), Tuple{ - DOF{Float64, Vertex}, # Temperature at vertices - DOF{Vec{3}, Vertex}, # Displacement at vertices - DOF{Float64, Cell}, # Pore pressure at cells - DOF{Float64, Edge} # Electric potential at edges - }} - - println("\nπŸ“‹ Field specification:") - println(" T: DOF{Float64, Vertex} - Temperature") - println(" u: DOF{Vec{3}, Vertex} - Displacement") - println(" p: DOF{Float64, Cell} - Pore pressure") - println(" Ο†: DOF{Float64, Edge} - Electric potential") - - # πŸš€ ONE ELEMENT CREATION CALL FOR ALL FIELDS! - println("\nπŸš€ Creating elements with ALL fields in ONE call...") - elements, mgr = create_elements!(mesh, Element{Tetrahedron{4}, Lagrange{1}, field_spec}) - - println(" βœ“ Created $(length(elements)) multi-field elements!") - println(" βœ“ Total DOFs: $(mgr.total_dofs)") - - # Extract DOF information from first element - elem = elements[1] - T_dofs = elem.dof_indices.T - u_dofs = elem.dof_indices.u - p_dofs = elem.dof_indices.p - Ο†_dofs = elem.dof_indices.Ο† - - println("\n✨ Element 1 DOF structure:") - println(" elem.dof_indices.T: $(T_dofs) ($(length(T_dofs)) DOFs)") - println(" elem.dof_indices.u: $(u_dofs) ($(length(u_dofs)) DOFs)") - println(" elem.dof_indices.p: $(p_dofs) ($(length(p_dofs)) DOFs)") - println(" elem.dof_indices.Ο†: $(Ο†_dofs) ($(length(Ο†_dofs)) DOFs)") - - # Verify DOF counts per element - @test length(T_dofs) == 4 # 4 vertices per tet - @test length(u_dofs) == 12 # 4 vertices Γ— 3 components - @test length(p_dofs) == 1 # 1 cell per element - @test length(Ο†_dofs) == 6 # 6 edges per tet - - println("\nπŸ“Š DOF Summary:") - println(" Total DOFs in system: $(mgr.total_dofs)") - - @test mgr.total_dofs == 26 # Verify total matches - - # Calculate DOF counts per field - n_T = 5 # 5 vertices (temperature) - n_u = 15 # 5 vertices Γ— 3 components (displacement) - n_p = 2 # 2 cells (pore pressure) - n_Ο† = 4 # 4 unique edges (electric potential) - n_total = n_T + n_u + n_p + n_Ο† # Should be 26 - - println("\n" * "="^70) - println("TOTAL SYSTEM:") - println("="^70) - println(" Temperature DOFs: $n_T") - println(" Displacement DOFs: $n_u (Vec{3})") - println(" Pore pressure DOFs: $n_p (Cell)") - println(" Electric DOFs: $n_Ο† (Edge)") - println(" " * "-"^40) - println(" TOTAL: $n_total") - - println("\nβœ“ Verification - Element DOF structure:") - println(" βœ… One element contains ALL four fields!") - println(" βœ… Type-safe field access: elem.dof_indices.T, .u, .p, .Ο†") - println(" βœ… T & u share vertex DOFs (natural coupling!)") - - @test length(elements[1].dof_indices.T) == 4 # 4 vertices - @test length(elements[1].dof_indices.u) == 12 # 4 vertices Γ— 3 components - @test length(elements[1].dof_indices.p) == 1 # 1 cell - @test length(elements[1].dof_indices.Ο†) == 6 # 6 edges - - println("\n" * "="^70) - println("ASSEMBLING COUPLED THM-E SYSTEM...") - println("="^70) - - # Material parameters - ΞΊ = 1.0 # Thermal conductivity - k = 1.0 # Hydraulic permeability - E = 1000.0 # Young's modulus - Ξ½ = 0.3 # Poisson's ratio - Ξ±_T = 1e-5 # Thermal expansion coefficient - Ξ±_p = 1e-3 # Poroelastic coefficient (Biot) - Οƒ_e = 1.0 # Electric conductivity - - println("\nπŸ“Š Material properties:") - println(" ΞΊ (thermal): $ΞΊ") - println(" k (hydraulic): $k") - println(" E (elastic): $E") - println(" Ξ½ (Poisson): $Ξ½") - println(" Ξ±_T (thermal): $Ξ±_T") - println(" Ξ±_p (Biot): $Ξ±_p") - println(" Οƒ_e (electric): $Οƒ_e") - - # Initialize block matrices for assembly - K_TT = spzeros(Float64, n_T, n_T) # Thermal diffusion - K_uu = spzeros(Float64, n_u, n_u) # Mechanical stiffness - K_uT = spzeros(Float64, n_u, n_T) # Thermal expansion coupling - K_up = spzeros(Float64, n_u, n_p) # Poroelastic coupling - K_pu = spzeros(Float64, n_p, n_u) # Consolidation coupling - K_pp = spzeros(Float64, n_p, n_p) # Hydraulic - K_φφ = spzeros(Float64, n_Ο†, n_Ο†) # Electric - - # Initialize global system matrix (full coupled system) - n_total = mgr.total_dofs - K_full = spzeros(Float64, n_total, n_total) - F_full = zeros(Float64, n_total) - - println("\nπŸ”§ Assembly strategy (COUPLED!):") - println(" 1. Thermal: K_TT from βˆ«ΞΊβˆ‡TΒ·βˆ‡T' dx") - println(" 2. Mechanical: K_uu from ∫C:Ξ΅(u):Ξ΅(u') dx") - println(" 3. Thermalβ†’Mech coupling: K_uT from ∫α_T*C:T*Ξ΅(u') dx") - println(" 4. Mechβ†’Pressure coupling: K_up from ∫α_p*p*βˆ‡Β·u' dx") - println(" 5. Pressureβ†’Mech: K_pu from βˆ«βˆ‡Β·u*p' dx (consolidation)") - println(" 6. Hydraulic: K_pp from ∫kβˆ‡pΒ·βˆ‡p' dx") - println(" 7. Electric: K_φφ from βˆ«Οƒ_eβˆ‡Ο†Β·βˆ‡Ο†' dx on edges") - println(" β†’ OFF-DIAGONAL blocks make this a TRULY COUPLED system!") - - # For simplicity: assemble diagonal blocks + key coupling terms - println("\nβš™οΈ Assembling COUPLED system (simplified for demo)...") - - # πŸš€ CRITICAL: Assembly loops iterate ONCE per element, accessing ALL fields! - # Each element contributes to MULTIPLE blocks simultaneously: - f_T = zeros(Float64, n_T) # Heat sources - f_u = zeros(Float64, n_u) # Body forces - f_p = zeros(Float64, n_p) # Fluid sources - f_Ο† = zeros(Float64, n_Ο†) # Charge sources - - # πŸš€ CRITICAL: Assembly loops iterate ONCE per element, accessing ALL fields! - # Each element contributes to MULTIPLE blocks simultaneously: - for elem in elements - # πŸŽ‰ ELEGANT: Extract DOFs for ALL fields from ONE element! - T_dofs_global = [Int(i) for i in elem.dof_indices.T] # Global DOF indices - u_dofs_global = [Int(i) for i in elem.dof_indices.u] - p_dof_global = Int(elem.dof_indices.p[1]) - Ο†_dofs_global = [Int(i) for i in elem.dof_indices.Ο†] - - # Map global DOFs to field-local indices (for block matrices) - # T field: DOFs 1-5 β†’ local 1-5 - T_dofs = T_dofs_global # Already 1-5 - - # u field: DOFs vary by node, but need to map to u-local indices - u_dofs_local = Int[] - for (i, g_dof) in enumerate(u_dofs_global) - # Find which local u DOF this is (1-based within u field) - # u field starts after T field - u_local = g_dof - n_T - if u_local > 0 && u_local <= n_u - push!(u_dofs_local, u_local) - end - end - - # p field: Cell DOF, need local index (1-2 for 2 cells) - p_dof_local = p_dof_global - (n_T + n_u) # Subtract T and u field sizes - - # Ο† field: Edge DOFs - Ο†_dofs_local = [g - (n_T + n_u + n_p) for g in Ο†_dofs_global] - - # 1. Thermal diffusion (diagonal) - for i in 1:4 - if T_dofs[i] > 0 && T_dofs[i] <= n_T - K_TT[T_dofs[i], T_dofs[i]] += ΞΊ * 0.1 - f_T[T_dofs[i]] += 0.01 # Heat source - end - end - - # 2. Mechanical stiffness (diagonal) - for u_local in u_dofs_local - if u_local > 0 && u_local <= n_u - K_uu[u_local, u_local] += E * 0.01 - end - end - - # 3. COUPLING: Thermal expansion (T β†’ u) - # K_uT couples displacement to temperature - for u_local in u_dofs_local, j in 1:length(T_dofs) - if u_local > 0 && u_local <= n_u && T_dofs[j] > 0 && T_dofs[j] <= n_T - K_uT[u_local, T_dofs[j]] += Ξ±_T * E * 0.001 # Mock coupling - end - end - - # 4. COUPLING: Poroelasticity (p β†’ u) - # K_up couples displacement to pressure - for u_local in u_dofs_local - if u_local > 0 && u_local <= n_u && p_dof_local > 0 && p_dof_local <= n_p - K_up[u_local, p_dof_local] += Ξ±_p * E * 0.002 # Mock coupling - end - end - - # 5. COUPLING: Consolidation (u β†’ p) - # K_pu couples pressure to displacement (symmetric) - for u_local in u_dofs_local - if p_dof_local > 0 && p_dof_local <= n_p && u_local > 0 && u_local <= n_u - K_pu[p_dof_local, u_local] += Ξ±_p * 0.002 # Mock coupling - end - end - - # 6. Pressure (diagonal) - if p_dof_local > 0 && p_dof_local <= n_p - K_pp[p_dof_local, p_dof_local] += k * 1.0 - end - - # 7. Electric (diagonal - edge basis) - for Ο†_local in Ο†_dofs_local - if Ο†_local > 0 && Ο†_local <= n_Ο† - K_φφ[Ο†_local, Ο†_local] += Οƒ_e * 0.05 - end - end - end - - println(" βœ“ All blocks assembled IN ONE PASS!") - println(" βœ“ Off-diagonal coupling terms included!") - println(" βœ“ This is TRUE multi-physics coupling!") - - # Build full coupled system - println("\nπŸ—οΈ Building COUPLED system matrix...") - - # DOF ranges for block assembly - T_dof_range = 1:n_T - u_dof_range = (n_T+1):(n_T+n_u) - p_dof_range = (n_T+n_u+1):(n_T+n_u+n_p) - Ο†_dof_range = (n_T+n_u+n_p+1):(n_T+n_u+n_p+n_Ο†) - - # Build block-by-block - K_full = spzeros(Float64, n_total, n_total) - - # Block (1,1): Thermal - K_full[T_dof_range, T_dof_range] = K_TT - - # Block (2,2): Mechanical - K_full[u_dof_range, u_dof_range] = K_uu - - # Block (2,1): Thermal-mechanical coupling - K_full[u_dof_range, T_dof_range] = K_uT - - # Block (2,3): Mechanical-pressure coupling - K_full[u_dof_range, p_dof_range] = K_up - - # Block (3,2): Pressure-mechanical coupling - K_full[p_dof_range, u_dof_range] = K_pu - - # Block (3,3): Pressure - K_full[p_dof_range, p_dof_range] = K_pp - - # Block (4,4): Electric - K_full[Ο†_dof_range, Ο†_dof_range] = K_φφ - - F_full = [f_T; f_u; f_p; f_Ο†] - - println(" System size: $(size(K_full))") - println(" Non-zeros: $(nnz(K_full))") - println(" Non-zeros in K_uT: $(nnz(K_uT)) ← Thermal-mechanical coupling!") - println(" Non-zeros in K_up: $(nnz(K_up)) ← Poroelastic coupling!") - println(" Non-zeros in K_pu: $(nnz(K_pu)) ← Consolidation coupling!") - println(" β†’ This is NOT block-diagonal! TRUE coupling!") - - # Verify coupling exists - @test nnz(K_uT) > 0 # Thermal expansion coupling must exist - # Note: K_up and K_pu may be zero in this simplified demo due to DOF layout - # @test nnz(K_up) > 0 # Poroelastic coupling must exist - # @test nnz(K_pu) > 0 # Consolidation coupling must exist - - # Apply boundary conditions - println("\nπŸ”’ Applying boundary conditions...") - # Fix enough DOFs to make system non-singular - # Fix all DOFs of first element to ensure solvability (this is a demo!) - bc_dofs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] # Fix T, u, and p for element 1 - for dof in bc_dofs - K_full[dof, :] .= 0.0 - K_full[:, dof] .= 0.0 - K_full[dof, dof] = 1.0 - F_full[dof] = 0.0 - end - println(" βœ“ Fixed $(length(bc_dofs)) DOFs (for demo purposes)") - - # Solve - println("\n🎯 SOLVING COUPLED THM-E SYSTEM...") - try - sol = K_full \ F_full - - # Extract fields (using correct total count) - T_sol = sol[1:5] # 5 T DOFs - u_p_Ο†_sol = sol[6:end] # Rest (u, p, Ο† mixed) - - println("\n" * "="^70) - println("✨ SOLUTION (showing first few DOFs):") - println("="^70) - - println("\nπŸ“Š Solution vector (first 10 DOFs):") - for i in 1:min(10, length(sol)) - println(" DOF $i: $(sol[i])") - end - - # Verify solution - @test all(isfinite.(sol)) - @test sol[1] β‰ˆ 0.0 atol=1e-10 # BC: T at node 1 = 0 - - println("\n βœ“ Solution obtained successfully!") - println(" βœ“ All values finite") - println(" βœ“ Boundary conditions satisfied") - - catch e - println("\n ⚠️ Solve failed (system may be under-constrained for full solve)") - println(" ⚠️ BUT: Assembly demonstrated successfully!") - @test true # Pass anyway - assembly is what matters - end - - println("\n" * "="^70) - println("πŸŽ‰ ACHIEVEMENTS UNLOCKED:") - println("="^70) - println(" βœ… Temperature DOFs at VERTICES") - println(" βœ… Displacement DOFs (Vec{3}) at VERTICES") - println(" βœ… Pore pressure DOFs at CELLS") - println(" βœ… Electric potential DOFs at EDGES") - println(" βœ… FOUR different entity types in ONE mesh!") - println(" βœ… GLOBAL DOF numbering across all fields") - println(" βœ… OFF-DIAGONAL coupling matrices (K_uT, K_up, K_pu)") - println(" βœ… SIMULTANEOUS assembly (one pass, all couplings!)") - println(" βœ… Full THM-E system solved ($n_total DOFs)") - println("\n πŸ† USE CASES:") - println(" β€’ Geothermal energy extraction") - println(" β€’ CO2 geological sequestration") - println(" β€’ Nuclear waste repository") - println(" β€’ Electrokinetic soil remediation") - println(" β€’ Hydraulic fracturing") - println(" β€’ Permafrost thawing") - println("="^70) - - println("\nπŸ’‘ THIS IS THE POWER OF THE CIARLET FRAMEWORK!") - println(" Different physics β†’ Different function spaces β†’ Different entities") - println(" BUT: All DOFs in ONE GLOBAL system β†’ TRUE coupling possible!") - println(" Assembly in ONE PASS β†’ Efficient and elegant! πŸš€") -end