replace Threads.nthreads by Threads.maxthreadid (#2612)

* replace Threads.nthreads by Threads.maxthreadid

* format

Author: ranocha

examples/tree_2d_dgsem/elixir_euler_vortex_amr.jl

@@ -15,7 +15,7 @@ function IndicatorVortex(semi)
     alpha = Vector{real(basis)}()
     A = Array{real(basis), 2}
     indicator_threaded = [A(undef, nnodes(basis), nnodes(basis))
-                          for _ in 1:Threads.nthreads()]
+                          for _ in 1:Threads.maxthreadid()]
     cache = (; semi.mesh, alpha, indicator_threaded)
 
     return IndicatorVortex{typeof(cache)}(cache)
@@ -119,7 +119,7 @@ initial_condition = initial_condition_isentropic_vortex
 # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
 # Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
 # To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
 # `StepsizeCallback` (CFL-Condition) and less diffusion.
 solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))
 
@@ -153,7 +153,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 # Add `:temperature` to `extra_node_variables` tuple ...
 extra_node_variables = (:temperature,)
 
-# ... and specify the function `get_node_variable` for this symbol, 
+# ... and specify the function `get_node_variable` for this symbol,
 # with first argument matching the symbol (turned into a type via `Val`) for dispatching.
 function Trixi.get_node_variable(::Val{:temperature}, u, mesh, equations, dg, cache)
     n_nodes = nnodes(dg)

src/solvers/dgmulti/dg.jl

@@ -189,7 +189,7 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
 
     # local storage for volume integral and source computations
     local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
+                             for _ in 1:Threads.maxthreadid()]
 
     # For curved meshes, we interpolate geometric terms from nodal points to quadrature points.
     # For affine meshes, we just access one element of this interpolated data.
@@ -200,9 +200,9 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
 
     # for scaling by curved geometric terms (not used by affine DGMultiMesh)
     flux_threaded = [[allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                      for _ in 1:NDIMS] for _ in 1:Threads.nthreads()]
+                      for _ in 1:NDIMS] for _ in 1:Threads.maxthreadid()]
     rotated_flux_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
+                             for _ in 1:Threads.maxthreadid()]
 
     return (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj,
             u_values, u_face_values, flux_face_values,

src/solvers/dgmulti/dg_parabolic.jl

@@ -38,13 +38,13 @@ function create_cache_parabolic(mesh::DGMultiMesh,
     gradients_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh))
 
     local_u_values_threaded = [similar(u_transformed, dg.basis.Nq)
-                               for _ in 1:Threads.nthreads()]
+                               for _ in 1:Threads.maxthreadid()]
     local_flux_viscous_threaded = [SVector{ndims(mesh)}(ntuple(_ -> similar(u_transformed,
                                                                             dg.basis.Nq),
                                                                ndims(mesh)))
-                                   for _ in 1:Threads.nthreads()]
+                                   for _ in 1:Threads.maxthreadid()]
     local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1])
-                                       for _ in 1:Threads.nthreads()]
+                                       for _ in 1:Threads.maxthreadid()]
 
     return (; u_transformed, gradients, flux_viscous,
             weak_differentiation_matrices, strong_differentiation_matrices,
@@ -152,9 +152,9 @@ function calc_gradient_interface_flux!(scalar_flux_face_values,
         idM, idP = mapM[face_node_index], mapP[face_node_index]
         uM = u_face_values[idM]
         uP = u_face_values[idP]
-        # Here, we use the "strong" formulation to compute the gradient. 
-        # This guarantees that the parabolic formulation is symmetric and 
-        # stable on curved meshes with variable geometric terms. 
+        # Here, we use the "strong" formulation to compute the gradient.
+        # This guarantees that the parabolic formulation is symmetric and
+        # stable on curved meshes with variable geometric terms.
         scalar_flux_face_values[idM] = 0.5f0 * (uP - uM)
     end
 

src/solvers/dgmulti/flux_differencing.jl

@@ -311,11 +311,11 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP, Real
     lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators
 
     local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
+                             for _ in 1:Threads.maxthreadid()]
 
     # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
     fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq)
-                               for _ in 1:Threads.nthreads()]
+                               for _ in 1:Threads.maxthreadid()]
 
     return (; md, Qrst_skew, dxidxhatj = md.rstxyzJ,
             invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
@@ -356,16 +356,16 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiff, RealT,
 
     # local storage for interface fluxes, rhs, and source
     local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
+                             for _ in 1:Threads.maxthreadid()]
 
     # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
     # The result is then transferred to rhs_local_threaded::StructArray{<:SVector} before
     # projecting it and storing it into `du`.
     fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, num_quad_points_total)
-                               for _ in 1:Threads.nthreads()]
+                               for _ in 1:Threads.maxthreadid()]
     rhs_local_threaded = [allocate_nested_array(uEltype, nvars,
                                                 (num_quad_points_total,), dg)
-                          for _ in 1:Threads.nthreads()]
+                          for _ in 1:Threads.maxthreadid()]
 
     # interpolate geometric terms to both quadrature and face values for curved meshes
     (; Vq, Vf) = dg.basis

src/solvers/dgmulti/flux_differencing_gauss_sbp.jl

@@ -427,9 +427,9 @@ function create_cache(mesh::DGMultiMesh, equations,
 
     nvars = nvariables(equations)
     rhs_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                                 for _ in 1:Threads.nthreads()]
+                                 for _ in 1:Threads.maxthreadid()]
     gauss_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                                   for _ in 1:Threads.nthreads()]
+                                   for _ in 1:Threads.maxthreadid()]
 
     return (; cache..., projection_matrix_gauss_to_face, gauss_LIFT, inv_gauss_weights,
             rhs_volume_local_threaded, gauss_volume_local_threaded,

src/solvers/dgmulti/shock_capturing.jl

@@ -39,8 +39,8 @@ function create_cache(::Type{IndicatorHennemannGassner}, equations::AbstractEqua
     alpha_tmp = similar(alpha)
 
     A = Vector{real(basis)}
-    indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]
-    modal_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]
+    indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
+    modal_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
 
     # initialize inverse Vandermonde matrices at Gauss-Legendre nodes
     (; N) = basis

src/solvers/dgmulti/types.jl

@@ -72,19 +72,19 @@ Base.real(rd::RefElemData) = eltype(rd.r)
 
 Create a discontinuous Galerkin method which uses
 - Approximations of polynomial degree `polydeg`.
-- Element type `element_type` (`Tri()`, `Quad()`, `Tet()`, `Hex()`, and `Wedge()` are 
+- Element type `element_type` (`Tri()`, `Quad()`, `Tet()`, `Hex()`, and `Wedge()` are
   currently supported)
 
 Optional:
-- `approximation_type` (default is `Polynomial()`; `SBP()` also supported for `Tri()`, 
+- `approximation_type` (default is `Polynomial()`; `SBP()` also supported for `Tri()`,
   `Quad()`, and `Hex()` element types).
-- `RefElemData_kwargs` are additional keyword arguments for `RefElemData`, such as 
+- `RefElemData_kwargs` are additional keyword arguments for `RefElemData`, such as
   `quad_rule_vol`.
-  
+
 For more info, see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/).
 
 !!! note "Wedge elements"
-    For `Wedge` elements (i.e. triangular prisms), the polynomial degree may optionally be 
+    For `Wedge` elements (i.e. triangular prisms), the polynomial degree may optionally be
     specified as a tuple of the form `polydeg = (polydeg_tri, polydeg_line)`.
 """
 function DGMulti(; polydeg = nothing,
@@ -100,9 +100,9 @@ function DGMulti(; polydeg = nothing,
             polydeg = polydeg, kwargs...)
 end
 
-# `Wedge` element types can optionally take `polydeg = (polydeg_tri, polydeg_line)`, which 
+# `Wedge` element types can optionally take `polydeg = (polydeg_tri, polydeg_line)`, which
 # constructs a `TensorProductWedge` approximation. Since Julia does not dispatch on keyword
-# arguments, we wrap a method which makes `polydeg` a positional argument. 
+# arguments, we wrap a method which makes `polydeg` a positional argument.
 function DGMulti(element_type::Wedge,
                  approximation_type,
                  volume_integral,
@@ -376,7 +376,7 @@ end
 function SimpleKronecker(NDIMS, A, eltype_A = eltype(A))
     @assert size(A, 1) == size(A, 2) # check if square
     tmp_storage = [zeros(eltype_A, ntuple(_ -> size(A, 2), NDIMS)...)
-                   for _ in 1:Threads.nthreads()]
+                   for _ in 1:Threads.maxthreadid()]
     return SimpleKronecker{NDIMS, typeof(A), typeof(tmp_storage)}(A, tmp_storage)
 end
 

src/solvers/dgsem_p4est/dg_2d.jl

@@ -14,11 +14,11 @@ function create_cache(mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}
     MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)},
                   uEltype, 2,
                   nvariables(equations) * nnodes(mortar_l2)}
-    fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
-    fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
-    fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
-    fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
-    u_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
+    fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
+    fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
+    fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
+    fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
+    u_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
 
     cache = (; fstar_primary_upper_threaded, fstar_primary_lower_threaded,
              fstar_secondary_upper_threaded, fstar_secondary_lower_threaded,

src/solvers/dgsem_p4est/dg_3d.jl

@@ -13,18 +13,18 @@ function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
     fstar_primary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
                                                 nnodes(mortar_l2),
                                                 nnodes(mortar_l2), 4)
-                              for _ in 1:Threads.nthreads()] |> VecOfArrays
+                              for _ in 1:Threads.maxthreadid()] |> VecOfArrays
     fstar_secondary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
                                                   nnodes(mortar_l2),
                                                   nnodes(mortar_l2), 4)
-                                for _ in 1:Threads.nthreads()] |> VecOfArrays
+                                for _ in 1:Threads.maxthreadid()] |> VecOfArrays
 
     fstar_tmp_threaded = [Array{uEltype, 3}(undef, nvariables(equations),
                                             nnodes(mortar_l2), nnodes(mortar_l2))
-                          for _ in 1:Threads.nthreads()] |> VecOfArrays
+                          for _ in 1:Threads.maxthreadid()] |> VecOfArrays
     u_threaded = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2),
                                     nnodes(mortar_l2))
-                  for _ in 1:Threads.nthreads()] |> VecOfArrays
+                  for _ in 1:Threads.maxthreadid()] |> VecOfArrays
 
     (; fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded, u_threaded)
 end

src/solvers/dgsem_tree/dg_1d.jl

@@ -42,9 +42,9 @@ function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations,
 
     A2dp1_x = Array{uEltype, 2}
     fstar1_L_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
-                                for _ in 1:Threads.nthreads()]
+                                for _ in 1:Threads.maxthreadid()]
     fstar1_R_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
-                                for _ in 1:Threads.nthreads()]
+                                for _ in 1:Threads.maxthreadid()]
 
     return (; cache..., fstar1_L_threaded, fstar1_R_threaded)
 end
@@ -54,9 +54,9 @@ function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations,
                       dg::DG, uEltype)
     A2dp1_x = Array{uEltype, 2}
     fstar1_L_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
-                                for _ in 1:Threads.nthreads()]
+                                for _ in 1:Threads.maxthreadid()]
     fstar1_R_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
-                                for _ in 1:Threads.nthreads()]
+                                for _ in 1:Threads.maxthreadid()]
 
     return (; fstar1_L_threaded, fstar1_R_threaded)
 end
@@ -347,7 +347,7 @@ end
     fstar1_R[:, nnodes(dg) + 1] .= zero(eltype(fstar1_R))
 
     for i in 2:nnodes(dg) # We compute FV02 fluxes at the (nnodes(dg) - 1) subcell boundaries
-        #             Reference element:             
+        #             Reference element:
         #  -1 ------------------0------------------ 1 -> x
         # Gauss-Lobatto-Legendre nodes (schematic for k = 3):
         #   .          .                  .         .
@@ -367,9 +367,9 @@ end
         # piecewise linear solution in both subcells next to the subcell interface.
         # Since these subcell boundaries are not aligned with the DG nodes,
         # on each neighboring subcell two linear solutions are reconstructed => 4 point stencil.
-        # For the outer interfaces the stencil shrinks since we do not consider values 
+        # For the outer interfaces the stencil shrinks since we do not consider values
         # outside the element (this is a volume integral).
-        # 
+        #
         # The left subcell node values are labelled `_ll` (left-left) and `_lr` (left-right), while
         # the right subcell node values are labelled `_rl` (right-left) and `_rr` (right-right).