Linear elasticity 1D

examples/tree_1d_dgsem/elixir_linearelasticity_convergence.jl

@@ -0,0 +1,52 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear elasticity equations
+
+rho = 3.0 # The material density rho can be any positive integer to ensure periodicity of the solution.
+c1_squared = 1.0 # Required to be one for the initial condition to stay periodic.
+c1 = sqrt(c1_squared)
+equations = LinearElasticityEquations1D(rho, c1_squared, c1)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
+
+coordinate_min = 0.0
+coordinate_max = 1.0
+
+mesh = TreeMesh(coordinate_min, coordinate_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+initial_condition = initial_condition_convergence_test
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (energy_kinetic,
+                                                                 energy_internal,
+                                                                 entropy))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 42.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

examples/tree_1d_dgsem/elixir_linearelasticity_impact.jl

@@ -0,0 +1,104 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear elasticity equations
+
+# Material parameters corresponding to steel
+rho = 7800.0 # kg/m³
+lambda = 9.3288e10
+mu = lambda
+equations = LinearElasticityEquations1D(rho = rho, mu = mu, lambda = lambda)
+
+basis = LobattoLegendreBasis(5)
+
+# Use subcell shock capturing technique to weaken oscillations at discontinuities
+alpha_max = 0.25 # This controls the amount of dissipation added (larger = more dissipation)
+indicator_variable = velocity # We limit here based on velocity
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = alpha_max,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = indicator_variable)
+
+volume_flux = flux_central
+surface_flux = flux_lax_friedrichs
+
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinate_min = -1.0
+coordinate_max = 1.0
+
+# E.g. hammer impact width of 20 cm
+impact_width() = 0.2
+# Refine the impact region to be able to represent the initial condition without oscillations
+refinement_patches = ((type = "box",
+                       coordinates_min = (-impact_width() / 2,),
+                       coordinates_max = (impact_width() / 2,)),)
+
+# Make sure to turn periodicity explicitly off as special boundary conditions are specified
+mesh = TreeMesh(coordinate_min, coordinate_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10_000,
+                refinement_patches = refinement_patches,
+                periodicity = false)
+
+function initial_condition_gaussian_impact(x, t, equations::LinearElasticityEquations1D)
+    amplitude = 1e6 # 1 MPa stress
+    edge_smoothing = 0.01 # Controls sharpness of edges (smaller = sharper)
+
+    # Smooth rectangle using tanh transitions
+    # Left edge: transition from 0 to 1 at x = -width/2
+    left_transition = 0.5 * (1 + tanh((x[1] + impact_width() / 2) / edge_smoothing))
+    # Right edge: transition from 1 to 0 at x = +width/2  
+    right_transition = 0.5 * (1 - tanh((x[1] - impact_width() / 2) / edge_smoothing))
+
+    # Combine both transitions to create smooth rectangle
+    sigma = amplitude * left_transition * right_transition
+
+    v = 0 # Initial displacement velocity is zero
+    return SVector(v, sigma)
+end
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+boundary_condition_dirichlet = BoundaryConditionDirichlet(initial_condition_gaussian_impact)
+
+boundary_conditions = (x_neg = boundary_condition_dirichlet,
+                       x_pos = boundary_condition_dirichlet)
+
+initial_condition = initial_condition_gaussian_impact
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1e-4) # Relatively short simulation time due to high wave speeds (~5990 m/s)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK54();
+            dt = 42.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -179,6 +179,7 @@ export AcousticPerturbationEquations2D,
        PolytropicEulerEquations2D,
        TrafficFlowLWREquations1D,
        MaxwellEquations1D,
+       LinearElasticityEquations1D,
        PassiveTracerEquations
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,

src/equations/equations.jl

@@ -724,4 +724,8 @@ include("traffic_flow_lwr_1d.jl")
 abstract type AbstractMaxwellEquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
 include("maxwell_1d.jl")
+
+abstract type AbstractLinearElasticityEquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("linear_elasticity_1d.jl")
 end # @muladd

src/equations/inviscid_burgers_1d.jl

@@ -67,9 +67,6 @@ Source terms used for convergence tests in combination with
     return SVector(du)
 end
 
-# Pre-defined source terms should be implemented as
-# function source_terms_WHATEVER(u, x, t, equations::InviscidBurgersEquation1D)
-
 # Calculate 1D flux for a single point
 @inline function flux(u, orientation::Integer, equation::InviscidBurgersEquation1D)
     return SVector(0.5f0 * u[1]^2)

src/equations/linear_elasticity_1d.jl

@@ -0,0 +1,143 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    LinearElasticityEquations1D(;rho::Real, lambda::Real, mu::Real)
+
+Linear elasticity equations in one space dimension. The equations are given by
+```math
+\partial_t
+\begin{pmatrix}
+    v_1 \\ \sigma_{11}
+\end{pmatrix}
++
+\partial_x
+\begin{pmatrix}
+    -\frac{1}{\rho} \sigma_{11} \\ -\frac{\rho}{c_1^2} v_1
+\end{pmatrix}
+=
+\begin{pmatrix}
+    0 \\ 0
+\end{pmatrix}.
+```
+The variables are the deformation velocity `v_1` and the stress `\sigma_{11}`.
+
+The parameters are the constant density of the material `\rho`
+and the Lamé parameters `\lambda` and `\mu`.
+From these, one can compute the dilatational wave speed as
+```math
+c_1^2= \frac{\lambda + 2 * \mu}{\rho}
+```
+In one dimension the linear elasticity equations reduce to a wave equation.
+
+For reference, see
+- Aleksey Sikstel (2020)
+  Analysis and numerical methods for coupled hyperbolic conservation laws
+  [DOI: 10.18154/RWTH-2020-07821](https://doi.org/10.18154/RWTH-2020-07821)
+"""
+struct LinearElasticityEquations1D{RealT <: Real} <:
+       AbstractLinearElasticityEquations{1, 2}
+    rho::RealT # Constant density of the material
+    c1::RealT # Dilatational wave speed
+    E::RealT # Young's modulus
+end
+
+function LinearElasticityEquations1D(; rho::Real, mu::Real, lambda::Real)
+    if !(rho > 0)
+        throw(ArgumentError("Density rho must be positive."))
+    end
+    if !(mu > 0)
+        throw(ArgumentError("Shear modulus mu (second Lamé parameter) must be positive."))
+    end
+
+    c1_squared = (lambda + 2 * mu) / rho
+
+    # Young's modulus
+    # See for reference equation (11-11) in
+    # https://cns.gatech.edu/~predrag/courses/PHYS-4421-04/lautrup/book/elastic.pdf
+    E = mu * (3 * lambda + 2 * mu) / (lambda + mu)
+    return LinearElasticityEquations1D(rho, sqrt(c1_squared), E)
+end
+
+function varnames(::typeof(cons2cons), ::LinearElasticityEquations1D)
+    ("v1", "sigma11")
+end
+function varnames(::typeof(cons2prim), ::LinearElasticityEquations1D)
+    ("v1", "sigma11")
+end
+
+"""
+    initial_condition_convergence_test(x, t, equations::LinearElasticityEquations1D)
+
+A smooth initial condition used for convergence tests.
+This requires that the material parameters `rho` is a positive integer
+and `c1` is equal to one.
+"""
+function initial_condition_convergence_test(x, t,
+                                            equations::LinearElasticityEquations1D)
+    @unpack rho = equations
+
+    v1 = sinpi(2 * t) * cospi(2 * x[1] / rho)
+    sigma11 = -cospi(2 * t) * sinpi(2 * x[1] * rho)
+
+    return SVector(v1, sigma11)
+end
+
+# Calculate 1D flux for a single point
+@inline function flux(u, orientation::Integer, equations::LinearElasticityEquations1D)
+    @unpack rho, c1 = equations
+    v1, sigma11 = u
+    f1 = -sigma11 / rho
+    f2 = -rho * c1^2 * v1
+
+    return SVector(f1, f2)
+end
+
+@inline have_constant_speed(::LinearElasticityEquations1D) = True()
+
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::LinearElasticityEquations1D)
+    return equations.c1
+end
+
+# Required for `StepsizeCallback`
+@inline function max_abs_speeds(equations::LinearElasticityEquations1D)
+    return equations.c1
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::LinearElasticityEquations1D)
+    @unpack c1 = equations
+
+    λ_min = -c1
+    λ_max = c1
+
+    return λ_min, λ_max
+end
+
+# Convert conservative variables to primitive
+@inline cons2prim(u, equations::LinearElasticityEquations1D) = u
+@inline cons2entropy(u, ::LinearElasticityEquations1D) = u
+
+# Useful for e.g. limiting indicator variable selection
+@inline function velocity(u, equations::LinearElasticityEquations1D)
+    return u[1]
+end
+
+@inline function energy_kinetic(u, equations::LinearElasticityEquations1D)
+    return 0.5f0 * equations.rho * u[1]^2
+end
+@inline function energy_internal(u, equations::LinearElasticityEquations1D)
+    return 0.5f0 * u[2]^2 / equations.E
+end
+@inline function energy_total(u, equations::LinearElasticityEquations1D)
+    return energy_kinetic(u, equations) + energy_internal(u, equations)
+end
+
+@inline entropy(u, equations::LinearElasticityEquations1D) = energy_total(u, equations)
+end # muladd

test/test_tree_1d.jl

@@ -52,6 +52,9 @@ isdir(outdir) && rm(outdir, recursive = true)
     # Maxwell
     include("test_tree_1d_maxwell.jl")
 
+    # Linear elasticity
+    include("test_tree_1d_linear_elasticity.jl")
+
     # Passive tracers
     include("test_tree_1d_passive_tracers.jl")
 end

test/test_tree_1d_linear_elasticity.jl

@@ -0,0 +1,40 @@
+module TestExamples1DLinearElasticity
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = joinpath(examples_dir(), "tree_1d_dgsem")
+
+@testset "Linear Elasticity" begin
+#! format: noindent
+
+@trixi_testset "elixir_linearelasticity_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_linearelasticity_convergence.jl"),
+                        analysis_callback=AnalysisCallback(semi,
+                                                           interval = analysis_interval,
+                                                           extra_analysis_errors = (:l2_error_primitive,
+                                                                                    :linf_error_primitive),
+                                                           extra_analysis_integrals = (energy_kinetic,
+                                                                                       energy_internal,
+                                                                                       entropy)),
+                        l2=[0.0007205516785218745, 0.0008036755866155103],
+                        linf=[0.0011507266875070855, 0.003249818227066381])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
+@trixi_testset "elixir_linearelasticity_impact.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearelasticity_impact.jl"),
+                        l2=[0.004334150310828556, 368790.1916121487],
+                        linf=[0.01070558926301203, 999999.9958777003])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+end
+
+end # module