Merge branch 'main' into remove_normal_direction_ll

trixi-framework · Sep 12, 2024 · c69719b · c69719b
2 parents 47dab1a + e4040e7
commit c69719b
Show file tree

Hide file tree

Showing 76 changed files with 1,039 additions and 897 deletions.
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
           julia  -e 'using JuliaFormatter; format(["benchmark", "examples", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |

diff --git a/Project.toml b/Project.toml
@@ -1,9 +1,10 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.8.9-DEV"
+version = "0.8.10-DEV"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -60,6 +61,7 @@ TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
 
 [compat]
+Accessors = "0.1.12"
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
 Convex = "0.16"

diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md
@@ -51,7 +51,7 @@ PRs that verify that running JuliaFormatter.jl again will not change the source
 To format your contributions before created a PR (or, at least, before requesting a review
 of your PR), you need to install JuliaFormatter.jl first by running
 ```shell
-julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
 ```
 You can then recursively format the core Julia files in the Trixi.jl repo by executing
 ```shell

diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -18,6 +18,7 @@ module Trixi
 # Include other packages that are used in Trixi.jl
 # (standard library packages first, other packages next, all of them sorted alphabetically)
 
+using Accessors: @reset
 using LinearAlgebra: LinearAlgebra, Diagonal, diag, dot, mul!, norm, cross, normalize, I,
                      UniformScaling, det
 using Printf: @printf, @sprintf, println

diff --git a/src/callbacks_step/alive.jl b/src/callbacks_step/alive.jl
@@ -45,7 +45,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         alive_callback = cb.affect!
 
         setup = [
-            "interval" => alive_callback.alive_interval,
+            "interval" => alive_callback.alive_interval
         ]
         summary_box(io, "AliveCallback", setup)
     end

diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
@@ -29,10 +29,39 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2},
     return (; u_local, u_tmp1, x_local, x_tmp1)
 end
 
+# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
+function create_cache_analysis(analyzer, mesh::P4estMesh{2, NDIMS_AMBIENT},
+                               equations, dg::DG, cache,
+                               RealT, uEltype) where {NDIMS_AMBIENT}
+
+    # pre-allocate buffers
+    # We use `StrideArray`s here since these buffers are used in performance-critical
+    # places and the additional information passed to the compiler makes them faster
+    # than native `Array`s.
+    u_local = StrideArray(undef, uEltype,
+                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)))
+    u_tmp1 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)))
+    x_local = StrideArray(undef, RealT,
+                          StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)))
+    x_tmp1 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)))
+    jacobian_local = StrideArray(undef, RealT,
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)))
+    jacobian_tmp1 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+
+    return (; u_local, u_tmp1, x_local, x_tmp1, jacobian_local, jacobian_tmp1)
+end
+
 function create_cache_analysis(analyzer,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
-                                           UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}},
+                                           UnstructuredMesh2D, T8codeMesh{2}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 

diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl
@@ -35,9 +35,51 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{3},
     return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2)
 end
 
+# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
 function create_cache_analysis(analyzer,
-                               mesh::Union{StructuredMesh{3}, P4estMesh{3},
-                                           T8codeMesh{3}},
+                               mesh::P4estMesh{3, NDIMS_AMBIENT},
+                               equations, dg::DG, cache,
+                               RealT, uEltype) where {NDIMS_AMBIENT}
+
+    # pre-allocate buffers
+    # We use `StrideArray`s here since these buffers are used in performance-critical
+    # places and the additional information passed to the compiler makes them faster
+    # than native `Array`s.
+    u_local = StrideArray(undef, uEltype,
+                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
+    u_tmp1 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
+    u_tmp2 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+    x_local = StrideArray(undef, RealT,
+                          StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
+    x_tmp1 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
+    x_tmp2 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+    jacobian_local = StrideArray(undef, RealT,
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)))
+    jacobian_tmp1 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)),
+                                StaticInt(nnodes(dg)))
+    jacobian_tmp2 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)),
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+
+    return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local,
+            jacobian_tmp1, jacobian_tmp2)
+end
+
+function create_cache_analysis(analyzer,
+                               mesh::Union{StructuredMesh{3}, T8codeMesh{3}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 

diff --git a/src/callbacks_step/averaging.jl b/src/callbacks_step/averaging.jl
@@ -45,7 +45,7 @@ function Base.show(io::IO, ::MIME"text/plain",
 
         setup = [
             "Start time" => first(averaging_callback.tspan),
-            "Final time" => last(averaging_callback.tspan),
+            "Final time" => last(averaging_callback.tspan)
         ]
         summary_box(io, "AveragingCallback", setup)
     end

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
@@ -48,7 +48,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-            "Selected semidiscretizations" => glm_speed_callback.semi_indices,
+            "Selected semidiscretizations" => glm_speed_callback.semi_indices
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end
@@ -83,7 +83,10 @@ function update_cleaning_speed!(semi, glm_speed_callback, dt)
     c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
 
     # c_h is proportional to its own time step divided by the complete MHD time step
-    equations.c_h = glm_scale * c_h_deltat / dt
+    # We use @reset here since the equations are immutable (to work on GPUs etc.).
+    # Thus, we need to modify the equations field of the semidiscretization.
+    @reset equations.c_h = glm_scale * c_h_deltat / dt
+    semi.equations = equations
 
     return semi
 end

diff --git a/src/callbacks_step/save_restart.jl b/src/callbacks_step/save_restart.jl
@@ -38,7 +38,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "interval" => save_restart_callback.interval,
             "save final solution" => save_restart_callback.save_final_restart ? "yes" :
                                      "no",
-            "output directory" => abspath(normpath(save_restart_callback.output_directory)),
+            "output directory" => abspath(normpath(save_restart_callback.output_directory))
         ]
         summary_box(io, "SaveRestartCallback", setup)
     end

diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl
@@ -62,7 +62,7 @@ function Base.show(io::IO, ::MIME"text/plain",
                                        "yes" : "no",
             "save final solution" => save_solution_callback.save_final_solution ?
                                      "yes" : "no",
-            "output directory" => abspath(normpath(save_solution_callback.output_directory)),
+            "output directory" => abspath(normpath(save_solution_callback.output_directory))
         ]
         summary_box(io, "SaveSolutionCallback", setup)
     end
@@ -85,7 +85,7 @@ function Base.show(io::IO, ::MIME"text/plain",
                                        "yes" : "no",
             "save final solution" => save_solution_callback.save_final_solution ?
                                      "yes" : "no",
-            "output directory" => abspath(normpath(save_solution_callback.output_directory)),
+            "output directory" => abspath(normpath(save_solution_callback.output_directory))
         ]
         summary_box(io, "SaveSolutionCallback", setup)
     end

diff --git a/src/callbacks_step/steady_state.jl b/src/callbacks_step/steady_state.jl
@@ -43,7 +43,7 @@ function Base.show(io::IO, ::MIME"text/plain",
 
         setup = [
             "absolute tolerance" => steady_state_callback.abstol,
-            "relative tolerance" => steady_state_callback.reltol,
+            "relative tolerance" => steady_state_callback.reltol
         ]
         summary_box(io, "SteadyStateCallback", setup)
     end

diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl
@@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         stepsize_callback = cb.affect!
 
         setup = [
-            "CFL number" => stepsize_callback.cfl_number,
+            "CFL number" => stepsize_callback.cfl_number
         ]
         summary_box(io, "StepsizeCallback", setup)
     end

diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl
@@ -72,7 +72,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "interval" => time_series_callback.interval,
             "solution_variables" => time_series_callback.solution_variables,
             "output_directory" => time_series_callback.output_directory,
-            "filename" => time_series_callback.filename,
+            "filename" => time_series_callback.filename
         ]
         summary_box(io, "TimeSeriesCallback", setup)
     end

diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl
@@ -50,7 +50,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "variable names" => visualization_callback.variable_names,
             "show mesh" => visualization_callback.show_mesh,
             "plot creator" => visualization_callback.plot_creator,
-            "plot data creator" => visualization_callback.plot_data_creator,
+            "plot data creator" => visualization_callback.plot_data_creator
         ]
         summary_box(io, "VisualizationCallback", setup)
     end

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
@@ -11,8 +11,8 @@
 The ideal compressible GLM-MHD equations for an ideal gas with ratio of
 specific heats `gamma` in two space dimensions.
 """
-mutable struct IdealGlmMhdEquations2D{RealT <: Real} <:
-               AbstractIdealGlmMhdEquations{2, 9}
+struct IdealGlmMhdEquations2D{RealT <: Real} <:
+       AbstractIdealGlmMhdEquations{2, 9}
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
     c_h::RealT                 # GLM cleaning speed
@@ -28,6 +28,11 @@ function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN)
     IdealGlmMhdEquations2D(promote(gamma, initial_c_h)...)
 end
 
+# Outer constructor for `@reset` works correctly
+function IdealGlmMhdEquations2D(gamma, inv_gamma_minus_one, c_h)
+    IdealGlmMhdEquations2D(gamma, c_h)
+end
+
 have_nonconservative_terms(::IdealGlmMhdEquations2D) = True()
 n_nonconservative_terms(::IdealGlmMhdEquations2D) = 2
 

diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
@@ -11,8 +11,8 @@
 The ideal compressible GLM-MHD equations for an ideal gas with ratio of
 specific heats `gamma` in three space dimensions.
 """
-mutable struct IdealGlmMhdEquations3D{RealT <: Real} <:
-               AbstractIdealGlmMhdEquations{3, 9}
+struct IdealGlmMhdEquations3D{RealT <: Real} <:
+       AbstractIdealGlmMhdEquations{3, 9}
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
     c_h::RealT                 # GLM cleaning speed
@@ -28,6 +28,11 @@ function IdealGlmMhdEquations3D(gamma; initial_c_h = convert(typeof(gamma), NaN)
     IdealGlmMhdEquations3D(promote(gamma, initial_c_h)...)
 end
 
+# Outer constructor for `@reset` works correctly
+function IdealGlmMhdEquations3D(gamma, inv_gamma_minus_one, c_h)
+    IdealGlmMhdEquations3D(gamma, c_h)
+end
+
 have_nonconservative_terms(::IdealGlmMhdEquations3D) = True()
 function varnames(::typeof(cons2cons), ::IdealGlmMhdEquations3D)
     ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3", "psi")

diff --git a/src/equations/ideal_glm_mhd_multicomponent_1d.jl b/src/equations/ideal_glm_mhd_multicomponent_1d.jl
@@ -10,8 +10,8 @@
 
 The ideal compressible multicomponent GLM-MHD equations in one space dimension.
 """
-mutable struct IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT <: Real} <:
-               AbstractIdealGlmMhdMulticomponentEquations{1, NVARS, NCOMP}
+struct IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT <: Real} <:
+       AbstractIdealGlmMhdMulticomponentEquations{1, NVARS, NCOMP}
     gammas::SVector{NCOMP, RealT}
     gas_constants::SVector{NCOMP, RealT}
     cv::SVector{NCOMP, RealT}
@@ -51,6 +51,11 @@ function IdealGlmMhdMulticomponentEquations1D(; gammas, gas_constants)
                                                                      __gas_constants)
 end
 
+# Outer constructor for `@reset` works correctly
+function IdealGlmMhdMulticomponentEquations1D(gammas, gas_constants, cv, cp, c_h)
+    IdealGlmMhdMulticomponentEquations1D(gammas = gammas, gas_constants = gas_constants)
+end
+
 @inline function Base.real(::IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT}) where {
                                                                                                NVARS,
                                                                                                NCOMP,

diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl
@@ -10,8 +10,8 @@
 
 The ideal compressible multicomponent GLM-MHD equations in two space dimensions.
 """
-mutable struct IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real} <:
-               AbstractIdealGlmMhdMulticomponentEquations{2, NVARS, NCOMP}
+struct IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real} <:
+       AbstractIdealGlmMhdMulticomponentEquations{2, NVARS, NCOMP}
     gammas::SVector{NCOMP, RealT}
     gas_constants::SVector{NCOMP, RealT}
     cv::SVector{NCOMP, RealT}
@@ -21,18 +21,19 @@ mutable struct IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real}
     function IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT}(gammas::SVector{NCOMP,
                                                                                        RealT},
                                                                        gas_constants::SVector{NCOMP,
-                                                                                              RealT}) where {
-                                                                                                             NVARS,
-                                                                                                             NCOMP,
-                                                                                                             RealT <:
-                                                                                                             Real
-                                                                                                             }
+                                                                                              RealT},
+                                                                       c_h::RealT) where {
+                                                                                          NVARS,
+                                                                                          NCOMP,
+                                                                                          RealT <:
+                                                                                          Real
+                                                                                          }
         NCOMP >= 1 ||
             throw(DimensionMismatch("`gammas` and `gas_constants` have to be filled with at least one value"))
 
         cv = gas_constants ./ (gammas .- 1)
         cp = gas_constants + gas_constants ./ (gammas .- 1)
-        c_h = convert(eltype(gammas), NaN)
+        c_h = convert(eltype(gammas), c_h)
 
         new(gammas, gas_constants, cv, cp, c_h)
     end
@@ -49,8 +50,30 @@ function IdealGlmMhdMulticomponentEquations2D(; gammas, gas_constants)
     __gammas = SVector(map(RealT, _gammas))
     __gas_constants = SVector(map(RealT, _gas_constants))
 
+    c_h = convert(RealT, NaN)
+
+    return IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT}(__gammas,
+                                                                     __gas_constants,
+                                                                     c_h)
+end
+
+# Outer constructor for `@reset` works correctly
+function IdealGlmMhdMulticomponentEquations2D(gammas, gas_constants, cv, cp, c_h)
+    _gammas = promote(gammas...)
+    _gas_constants = promote(gas_constants...)
+    RealT = promote_type(eltype(_gammas), eltype(_gas_constants))
+
+    NVARS = length(_gammas) + 8
+    NCOMP = length(_gammas)
+
+    __gammas = SVector(map(RealT, _gammas))
+    __gas_constants = SVector(map(RealT, _gas_constants))
+
+    c_h = convert(RealT, c_h)
+
     return IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT}(__gammas,
-                                                                     __gas_constants)
+                                                                     __gas_constants,
+                                                                     c_h)
 end
 
 @inline function Base.real(::IdealGlmMhdMulticomponentEquations2D{NVARS, NCOMP, RealT}) where {