Merge branch 'main' into ranocha-patch-2

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"

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

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)
 

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)
 

src/callbacks_step/glm_speed.jl

@@ -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

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
 

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")

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,

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 {

src/meshes/p4est_mesh.jl

@@ -6,12 +6,23 @@
 #! format: noindent
 
 """
-    P4estMesh{NDIMS} <: AbstractMesh{NDIMS}
+    P4estMesh{NDIMS, NDIMS_AMBIENT} <: AbstractMesh{NDIMS}
 
 An unstructured curved mesh based on trees that uses the C library `p4est`
 to manage trees and mesh refinement.
+
+The parameter `NDIMS` denotes the dimension of the spatial domain or manifold represented
+by the mesh itself, while `NDIMS_AMBIENT` denotes the dimension of the ambient space in
+which the mesh is embedded. For example, the type `P4estMesh{3, 3}` corresponds to a
+standard mesh for a three-dimensional volume, whereas `P4estMesh{2, 3}` corresponds to a
+mesh for a two-dimensional surface or shell in three-dimensional space.
+
+!!! warning "Experimental implementation"
+    The use of `NDIMS != NDIMS_AMBIENT` is an experimental feature and may change in future
+    releases.
 """
-mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <:
+mutable struct P4estMesh{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost,
+                         NDIMSP2, NNODES} <:
                AbstractMesh{NDIMS}
     p4est       :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t}
     is_parallel :: IsParallel
@@ -48,7 +59,14 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
         ghost = ghost_new_p4est(p4est)
         ghost_pw = PointerWrapper(ghost)
 
-        mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel),
+        # To enable the treatment of a manifold of dimension NDIMS embedded within an
+        # ambient space of dimension NDIMS_AMBIENT, we store both as type parameters and
+        # allow them to differ in the general case. This functionality is used for
+        # constructing discretizations on spherical shell domains for applications in
+        # global atmospheric modelling. The ambient dimension NDIMS_AMBIENT is therefore 
+        # set here in the inner constructor to size(tree_node_coordinates, 1).
+        mesh = new{NDIMS, size(tree_node_coordinates, 1),
+                   eltype(tree_node_coordinates), typeof(is_parallel),
                    typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw,
                                                                                  is_parallel,
                                                                                  ghost_pw,
@@ -66,8 +84,8 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
     end
 end
 
-const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:False}
-const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True}
+const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:False}
+const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:True}
 
 @inline mpi_parallel(mesh::SerialP4estMesh) = False()
 @inline mpi_parallel(mesh::ParallelP4estMesh) = True()
@@ -87,7 +105,8 @@ function destroy_mesh(mesh::P4estMesh{3})
 end
 
 @inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS
-@inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
+@inline Base.real(::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} = RealT
+@inline ndims_ambient(::P4estMesh{NDIMS, NDIMS_AMBIENT}) where {NDIMS, NDIMS_AMBIENT} = NDIMS_AMBIENT
 
 @inline function ntrees(mesh::P4estMesh)
     return mesh.p4est.trees.elem_count[]
@@ -97,7 +116,8 @@ end
 @inline ncellsglobal(mesh::P4estMesh) = Int(mesh.p4est.global_num_quadrants[])
 
 function Base.show(io::IO, mesh::P4estMesh)
-    print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}")
+    print(io, "P4estMesh{", ndims(mesh), ", ", ndims_ambient(mesh), ", ", real(mesh),
+          "}")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
@@ -110,8 +130,9 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
             "polydeg" => length(mesh.nodes) - 1
         ]
         summary_box(io,
-                    "P4estMesh{" * string(ndims(mesh)) * ", " * string(real(mesh)) *
-                    "}", setup)
+                    "P4estMesh{" * string(ndims(mesh)) * ", " *
+                    string(ndims_ambient(mesh)) *
+                    ", " * string(real(mesh)) * "}", setup)
     end
 end
 

src/semidiscretization/semidiscretization_coupled.jl

@@ -16,7 +16,8 @@ The semidiscretizations can be coupled by gluing meshes together using [`Boundar
 !!! warning "Experimental code"
     This is an experimental feature and can change any time.
 """
-struct SemidiscretizationCoupled{S, Indices, EquationList} <: AbstractSemidiscretization
+mutable struct SemidiscretizationCoupled{S, Indices, EquationList} <:
+               AbstractSemidiscretization
     semis::S
     u_indices::Indices # u_ode[u_indices[i]] is the part of u_ode corresponding to semis[i]
     performance_counter::PerformanceCounter
@@ -383,7 +384,10 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
         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
     end
 
     return semi_coupled