Replace mutable structs using Accessors

Project.toml

@@ -4,6 +4,7 @@ authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor
 version = "0.8.9-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"

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

src/semidiscretization/semidiscretization_hyperbolic.jl

@@ -11,10 +11,10 @@
 A struct containing everything needed to describe a spatial semidiscretization
 of a hyperbolic conservation law.
 """
-struct SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition,
-                                    BoundaryConditions,
-                                    SourceTerms, Solver, Cache} <:
-       AbstractSemidiscretization
+mutable struct SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition,
+                                            BoundaryConditions,
+                                            SourceTerms, Solver, Cache} <:
+               AbstractSemidiscretization
     mesh::Mesh
     equations::Equations