diff --git a/examples/tree_1d_dgsem/elixir_linearelasticity_convergence.jl b/examples/tree_1d_dgsem/elixir_linearelasticity_convergence.jl new file mode 100644 index 00000000000..0c7795c7af0 --- /dev/null +++ b/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); diff --git a/examples/tree_1d_dgsem/elixir_linearelasticity_impact.jl b/examples/tree_1d_dgsem/elixir_linearelasticity_impact.jl new file mode 100644 index 00000000000..4a875f0fab9 --- /dev/null +++ b/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); diff --git a/src/Trixi.jl b/src/Trixi.jl index 6b596bd74fb..a85ad5363e4 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -179,6 +179,7 @@ export AcousticPerturbationEquations2D, PolytropicEulerEquations2D, TrafficFlowLWREquations1D, MaxwellEquations1D, + LinearElasticityEquations1D, PassiveTracerEquations export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 850aec70ad5..c90549fd529 100644 --- a/src/equations/equations.jl +++ b/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 diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 8374ce04562..206ac8a94c0 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/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) diff --git a/src/equations/linear_elasticity_1d.jl b/src/equations/linear_elasticity_1d.jl new file mode 100644 index 00000000000..a7d2e782ee3 --- /dev/null +++ b/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 diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index e16d4bcd1bc..3be469b1bb4 100644 --- a/test/test_tree_1d.jl +++ b/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 diff --git a/test/test_tree_1d_linear_elasticity.jl b/test/test_tree_1d_linear_elasticity.jl new file mode 100644 index 00000000000..3d028d5e7a2 --- /dev/null +++ b/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