Transient multiphysics (time integration)

Overview

Parabolic and hyperbolic problems are marched in time with the native θ-method (transient_solve) and Newmark-β (newmark_solve). Both operate on the mass matrix M = ∫ uv and the stiffness K, assembled through the weak-form layer.

using Rasmah

mesh = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
U = FESpace(mesh, ReferenceFE(3, 1))

M = assemble_mass(U)
K = assemble_matrix(heat_form(U, thermal_conductivity(steel)))

125×125 SparseMatrixCSC{Float64}:
 12.975  -4.325    0.0      0.0     …   0.0      0.0      0.0      0.0
 -4.325  21.625   -4.325    0.0         0.0      0.0      0.0      0.0
  0.0    -4.325   21.625   -4.325       0.0      0.0      0.0      0.0
  0.0     0.0     -4.325   21.625       0.0      0.0      0.0      0.0
  0.0     0.0      0.0     -4.325       0.0      0.0      0.0      0.0
 -4.325   0.0      0.0      0.0     …   0.0      0.0      0.0      0.0
  0.0    -6.4875   0.0      0.0         0.0      0.0      0.0      0.0
  0.0     0.0     -6.4875   0.0         0.0      0.0      0.0      0.0
  0.0     0.0      0.0     -6.4875      0.0      0.0      0.0      0.0
  0.0     0.0      0.0      0.0         0.0      0.0      0.0      0.0
  ⋮                                 ⋱                             
  0.0     0.0      0.0      0.0        -6.4875   0.0      0.0      0.0
  0.0     0.0      0.0      0.0         0.0     -6.4875   0.0      0.0
  0.0     0.0      0.0      0.0         0.0      0.0     -6.4875   0.0
  0.0     0.0      0.0      0.0         0.0      0.0      0.0     -4.325
  0.0     0.0      0.0      0.0     …  -4.325    0.0      0.0      0.0
  0.0     0.0      0.0      0.0        21.625   -4.325    0.0      0.0
  0.0     0.0      0.0      0.0        -4.325   21.625   -4.325    0.0
  0.0     0.0      0.0      0.0         0.0     -4.325   21.625   -4.325
  0.0     0.0      0.0      0.0         0.0      0.0     -4.325   12.975

Transient heat (θ-method)

Cool a hot block with u = 0 on the x = 0 and x = 1 faces:

fixed = Int[]
for i in 1:U.nnodes
    x = U.nodes[:, i]
    (x[1] < 1e-9 || x[1] > 1 - 1e-9) && push!(fixed, i)
end
free = setdiff(1:ndofs(U), fixed)

u0 = ones(length(free))                                # hot interior
T = transient_solve(M[free, free], K[free, free], u0,
                    t -> zeros(length(free)), 100, 0.01)   # backward Euler

println("initial max T = ", maximum(T[1]), ", final max T = ", round(maximum(T[end]); digits = 3))

initial max T = 1.0, final max T = 0.0

θ = 1 is backward Euler (default), θ = 1/2 Crank–Nicolson, θ = 0 forward Euler.

Structural dynamics (Newmark-β)

M ü + K u = f(t) with the average-acceleration scheme (β = 1/4, γ = 1/2):

u0d = zeros(length(free))
v0  = zeros(length(free))
Uhist = newmark_solve(M[free, free], K[free, free], u0d, v0,
                      t -> zeros(length(free)), 50, 0.01)
println("Newmark steps: ", length(Uhist))

Newmark steps: 51

Both solvers are built on NativeLinearAlgebra factorizations, so they reuse the same preconditioner/direct machinery as the steady-state physics.