Merge branch 'trixi-framework:main' into structprint

examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl

@@ -136,7 +136,9 @@ function Trixi.analyze(::Val{:energy_potential}, du, u_euler, t,
         u_gravity_local = Trixi.get_node_vars(u_gravity, equations_gravity, dg, i, j,
                                               element)
         # OBS! subtraction is specific to Jeans instability test where rho0 = 1.5e7
-        return (u_euler_local[1] - 1.5e7) * u_gravity_local[1]
+        # For formula of potential energy see
+        # "Galactic Dynamics" by Binney and Tremaine, 2nd ed., equation (2.18)
+        return 0.5f0 * (u_euler_local[1] - 1.5f7) * u_gravity_local[1]
     end
     return e_potential
 end

src/equations/lattice_boltzmann_2d.jl

@@ -235,9 +235,6 @@ No-slip wall boundary condition using the bounce-back approach.
     return flux
 end
 
-# Pre-defined source terms should be implemented as
-# function source_terms_WHATEVER(u, x, t, equations::LatticeBoltzmannEquations2D)
-
 # Calculate 1D flux in for a single point
 @inline function flux(u, orientation::Integer, equations::LatticeBoltzmannEquations2D)
     if orientation == 1
@@ -248,11 +245,6 @@ end
     return v_alpha .* u
 end
 
-# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
-# @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::LatticeBoltzmannEquations2D)
-#   λ_max =
-# end
-
 """
     flux_godunov(u_ll, u_rr, orientation, 
                  equations::LatticeBoltzmannEquations2D)

src/equations/lattice_boltzmann_3d.jl

@@ -223,9 +223,6 @@ function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations3D
     return equilibrium_distribution(rho, v1, v2, v3, equations)
 end
 
-# Pre-defined source terms should be implemented as
-# function source_terms_WHATEVER(u, x, t, equations::LatticeBoltzmannEquations3D)
-
 # Calculate 1D flux in for a single point
 @inline function flux(u, orientation::Integer, equations::LatticeBoltzmannEquations3D)
     if orientation == 1 # x-direction
@@ -238,11 +235,6 @@ end
     return v_alpha .* u
 end
 
-# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
-# @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::LatticeBoltzmannEquations3D)
-#   λ_max =
-# end
-
 """
     flux_godunov(u_ll, u_rr, orientation, 
                  equations::LatticeBoltzmannEquations3D)

src/semidiscretization/semidiscretization_euler_gravity.jl

@@ -284,11 +284,6 @@ end
 function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
     @unpack semi_euler, semi_gravity, parameters, gravity_counter, cache = semi
 
-    # Can be changed by AMR
-    resize!(cache.du_ode, length(cache.u_ode))
-    resize!(cache.u_tmp1_ode, length(cache.u_ode))
-    resize!(cache.u_tmp2_ode, length(cache.u_ode))
-
     u_euler = wrap_array(u_ode, semi_euler)
     u_gravity = wrap_array(cache.u_ode, semi_gravity)
     du_gravity = wrap_array(cache.du_ode, semi_gravity)
@@ -568,7 +563,17 @@ end
                                              t, iter; kwargs...)
     passive_args = ((semi.cache.u_ode,
                      mesh_equations_solver_cache(semi.semi_gravity)...),)
-    amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)..., semi, t, iter;
-                 kwargs..., passive_args = passive_args)
+    has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)...,
+                               semi, t, iter;
+                               kwargs..., passive_args = passive_args)
+
+    if has_changed
+        new_length = length(semi.cache.u_ode)
+        resize!(semi.cache.du_ode, new_length)
+        resize!(semi.cache.u_tmp1_ode, new_length)
+        resize!(semi.cache.u_tmp2_ode, new_length)
+    end
+
+    return has_changed
 end
 end # @muladd

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -40,31 +40,27 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)
     nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT)
     inverse_weights_ = inv.(weights_)
 
-    _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_, RealT)
+    _, inverse_vandermonde_legendre = vandermonde_legendre(nodes_, RealT)
 
-    boundary_interpolation_ = zeros(RealT, nnodes_, 2)
-    boundary_interpolation_[:, 1] = calc_lhat(-one(RealT), nodes_, weights_)
-    boundary_interpolation_[:, 2] = calc_lhat(one(RealT), nodes_, weights_)
+    boundary_interpolation = zeros(RealT, nnodes_, 2)
+    boundary_interpolation[:, 1] = calc_lhat(-one(RealT), nodes_, weights_)
+    boundary_interpolation[:, 2] = calc_lhat(one(RealT), nodes_, weights_)
 
-    derivative_matrix_ = polynomial_derivative_matrix(nodes_)
-    derivative_split_ = calc_dsplit(nodes_, weights_)
-    derivative_split_transpose_ = Matrix(derivative_split_')
-    derivative_dhat_ = calc_dhat(nodes_, weights_)
+    derivative_matrix = polynomial_derivative_matrix(nodes_)
+    derivative_split = calc_dsplit(nodes_, weights_)
+    derivative_split_transpose = Matrix(derivative_split')
+    derivative_dhat = calc_dhat(nodes_, weights_)
 
-    # type conversions to get the requested real type and enable possible
-    # optimizations of runtime performance and latency
+    # Type conversions to enable possible optimizations of runtime performance 
+    # and latency
     nodes = SVector{nnodes_, RealT}(nodes_)
     weights = SVector{nnodes_, RealT}(weights_)
     inverse_weights = SVector{nnodes_, RealT}(inverse_weights_)
 
-    inverse_vandermonde_legendre = convert.(RealT, inverse_vandermonde_legendre_)
-    boundary_interpolation = convert.(RealT, boundary_interpolation_)
-
-    # Usually as fast as `SMatrix` (when using `let` in the volume integral/`@threaded`)
-    derivative_matrix = Matrix{RealT}(derivative_matrix_)
-    derivative_split = Matrix{RealT}(derivative_split_)
-    derivative_split_transpose = Matrix{RealT}(derivative_split_transpose_)
-    derivative_dhat = Matrix{RealT}(derivative_dhat_)
+    # We keep the matrices above stored using the standard `Matrix` type 
+    # since this is usually as fast as `SMatrix`
+    # (when using `let` in the volume integral/`@threaded`)
+    # and reduces latency
 
     return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes),
                                 typeof(inverse_vandermonde_legendre),
@@ -163,17 +159,15 @@ function MortarL2(basis::LobattoLegendreBasis)
     RealT = real(basis)
     nnodes_ = nnodes(basis)
 
-    forward_upper_ = calc_forward_upper(nnodes_, RealT)
-    forward_lower_ = calc_forward_lower(nnodes_, RealT)
-    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
-    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
-
-    # type conversions to get the requested real type and enable possible
-    # optimizations of runtime performance and latency
+    forward_upper = calc_forward_upper(nnodes_, RealT)
+    forward_lower = calc_forward_lower(nnodes_, RealT)
+    reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
+    reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
 
-    # Usually as fast as `SMatrix` but better for latency
-    forward_upper = Matrix{RealT}(forward_upper_)
-    forward_lower = Matrix{RealT}(forward_lower_)
+    # We keep the matrices above stored using the standard `Matrix` type 
+    # since this is usually as fast as `SMatrix`
+    # (when using `let` in the volume integral/`@threaded`)
+    # and reduces latency
 
     # TODO: Taal performance
     #       Check the performance of different implementations of `mortar_fluxes_to_elements!`
@@ -183,8 +177,6 @@ function MortarL2(basis::LobattoLegendreBasis)
     #       `@tullio` when the matrix sizes are not necessarily static.
     # reverse_upper = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_upper_)
     # reverse_lower = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_lower_)
-    reverse_upper = Matrix{RealT}(reverse_upper_)
-    reverse_lower = Matrix{RealT}(reverse_lower_)
 
     LobattoLegendreMortarL2{RealT, nnodes_, typeof(forward_upper),
                             typeof(reverse_upper)}(forward_upper, forward_lower,

src/solvers/dgsem/dgsem.jl

@@ -22,7 +22,7 @@ Create a discontinuous Galerkin spectral element method (DGSEM) using a
 """
 const DGSEM = DG{Basis} where {Basis <: LobattoLegendreBasis}
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(basis::LobattoLegendreBasis,
                surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm(),
@@ -32,7 +32,7 @@ function DGSEM(basis::LobattoLegendreBasis,
               typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral)
 end
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(basis::LobattoLegendreBasis,
                surface_integral::AbstractSurfaceIntegral,
                volume_integral = VolumeIntegralWeakForm(),
@@ -41,7 +41,7 @@ function DGSEM(basis::LobattoLegendreBasis,
               typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral)
 end
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(RealT, polydeg::Integer,
                surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm(),
@@ -51,6 +51,7 @@ function DGSEM(RealT, polydeg::Integer,
     return DGSEM(basis, surface_flux, volume_integral, mortar)
 end
 
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(polydeg, surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm())
     DGSEM(Float64, polydeg, surface_flux, volume_integral)

src/solvers/dgsem/l2projection.jl

@@ -29,9 +29,9 @@ function calc_forward_upper(n_nodes, RealT = Float64)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] + 1), nodes, wbary)
+        poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] + 1), nodes, wbary)
         for i in 1:n_nodes
             operator[j, i] = poly[i]
         end
@@ -49,9 +49,9 @@ function calc_forward_lower(n_nodes, RealT = Float64)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] - 1), nodes, wbary)
+        poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] - 1), nodes, wbary)
         for i in 1:n_nodes
             operator[j, i] = poly[i]
         end
@@ -70,12 +70,12 @@ function calc_reverse_upper(n_nodes, ::Val{:gauss}, RealT = Float64)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] + 1),
+        poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] + 1),
                                                   gauss_nodes, gauss_wbary)
         for i in 1:n_nodes
-            operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i]
+            operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i]
         end
     end
 
@@ -97,12 +97,12 @@ function calc_reverse_lower(n_nodes, ::Val{:gauss}, RealT = Float64)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] - 1),
+        poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] - 1),
                                                   gauss_nodes, gauss_wbary)
         for i in 1:n_nodes
-            operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i]
+            operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i]
         end
     end
 

src/solvers/dgsem_tree/dg_2d.jl

@@ -1118,7 +1118,7 @@ end
     #   @views mul!(surface_flux_values[v, :, direction, large_element],
     #               mortar_l2.reverse_upper, fstar_upper[v, :])
     #   @views mul!(surface_flux_values[v, :, direction, large_element],
-    #               mortar_l2.reverse_lower,  fstar_lower[v, :], true, true)
+    #               mortar_l2.reverse_lower, fstar_lower[v, :], true, true)
     # end
     # The code above could be replaced by the following code. However, the relative efficiency
     # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper.

src/solvers/dgsem_tree/dg_2d_parabolic.jl

@@ -725,7 +725,7 @@ end
     #   @views mul!(surface_flux_values[v, :, direction, large_element],
     #               mortar_l2.reverse_upper, fstar_upper[v, :])
     #   @views mul!(surface_flux_values[v, :, direction, large_element],
-    #               mortar_l2.reverse_lower,  fstar_lower[v, :], true, true)
+    #               mortar_l2.reverse_lower, fstar_lower[v, :], true, true)
     # end
     # The code above could be replaced by the following code. However, the relative efficiency
     # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper.