diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 8a11429766a..be9871666a8 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -710,18 +710,6 @@ function Base.show(io::IO, mime::MIME"text/plain", end end -# Required to be able to run `SimpleSSPRK33` without `VolumeIntegralSubcellLimiting` -Base.resize!(semi, volume_integral::AbstractVolumeIntegral, new_size) = nothing - -function Base.resize!(semi, volume_integral::VolumeIntegralSubcellLimiting, new_size) - # Resize container antidiffusive_fluxes - resize!(semi.cache.antidiffusive_fluxes, new_size) - - # Resize container subcell_limiter_coefficients - @unpack limiter = volume_integral - return resize!(limiter.cache.subcell_limiter_coefficients, new_size) -end - # TODO: FD. Should this definition live in a different file because it is # not strictly a DG method? """ diff --git a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl index ed57b91ce1c..00462d05786 100644 --- a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl @@ -24,7 +24,8 @@ function create_cache(mesh::P4estMesh{3}, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.maxthreadid()] - antidiffusive_fluxes = ContainerAntidiffusiveFlux3D{uEltype}(0, + n_elements = nelements(cache_containers.elements) + antidiffusive_fluxes = ContainerAntidiffusiveFlux3D{uEltype}(n_elements, nvariables(equations), nnodes(dg)) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index 4b1e10da09a..feb4a147136 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -537,7 +537,8 @@ end variable, min_or_max, initial_check, final_check, inverse_jacobian, dt, - equations, dg, cache, limiter) + equations::AbstractEquations{3}, + dg, cache, limiter) (; inverse_weights) = dg.basis # Plays role of inverse DG-subcell sizes (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index c7d7b71b432..6c5d970f70c 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -1342,6 +1342,9 @@ function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables, antidiffusive_flux2_R = unsafe_wrap(Array, pointer(_antidiffusive_flux2_R), (n_variables, n_nodes, n_nodes + 1, capacity)) + reset_antidiffusive_fluxes!(antidiffusive_flux1_L, antidiffusive_flux1_R, + antidiffusive_flux2_L, antidiffusive_flux2_R) + return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, @@ -1380,17 +1383,22 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity) (n_variables, n_nodes, n_nodes + 1, capacity)) - uEltype = eltype(fluxes.antidiffusive_flux1_L) - @threaded for element in axes(fluxes.antidiffusive_flux1_L, 4) - fluxes.antidiffusive_flux1_L[:, 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_L[:, n_nodes + 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_R[:, 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_R[:, n_nodes + 1, :, element] .= zero(uEltype) - - fluxes.antidiffusive_flux2_L[:, :, 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_L[:, :, n_nodes + 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_R[:, :, 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_R[:, :, n_nodes + 1, element] .= zero(uEltype) + return nothing +end + +function reset_antidiffusive_fluxes!(antidiffusive_flux1_L, antidiffusive_flux1_R, + antidiffusive_flux2_L, antidiffusive_flux2_R) + uEltype = eltype(antidiffusive_flux1_L) + @threaded for element in axes(antidiffusive_flux1_L, 4) + antidiffusive_flux1_L[:, 1, :, element] .= zero(uEltype) + antidiffusive_flux1_L[:, end, :, element] .= zero(uEltype) + antidiffusive_flux1_R[:, 1, :, element] .= zero(uEltype) + antidiffusive_flux1_R[:, end, :, element] .= zero(uEltype) + + antidiffusive_flux2_L[:, :, 1, element] .= zero(uEltype) + antidiffusive_flux2_L[:, :, end, element] .= zero(uEltype) + antidiffusive_flux2_R[:, :, 1, element] .= zero(uEltype) + antidiffusive_flux2_R[:, :, end, element] .= zero(uEltype) end return nothing diff --git a/src/solvers/dgsem_tree/containers_3d.jl b/src/solvers/dgsem_tree/containers_3d.jl index da398ef1938..09aebe049b5 100644 --- a/src/solvers/dgsem_tree/containers_3d.jl +++ b/src/solvers/dgsem_tree/containers_3d.jl @@ -845,6 +845,11 @@ function ContainerAntidiffusiveFlux3D{uEltype}(capacity::Integer, n_variables, antidiffusive_flux3_R = unsafe_wrap(Array, pointer(_antidiffusive_flux3_R), (n_variables, n_nodes, n_nodes, n_nodes + 1, capacity)) + + reset_antidiffusive_fluxes!(antidiffusive_flux1_L, antidiffusive_flux1_R, + antidiffusive_flux2_L, antidiffusive_flux2_R, + antidiffusive_flux3_L, antidiffusive_flux3_R) + return ContainerAntidiffusiveFlux3D{uEltype}(antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, @@ -908,22 +913,28 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux3D, capacity) n_nodes, n_nodes, n_nodes + 1, capacity)) - uEltype = eltype(fluxes.antidiffusive_flux1_L) - @threaded for element in axes(fluxes.antidiffusive_flux1_L, 5) - fluxes.antidiffusive_flux1_L[:, 1, :, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_L[:, n_nodes + 1, :, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_R[:, 1, :, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux1_R[:, n_nodes + 1, :, :, element] .= zero(uEltype) - - fluxes.antidiffusive_flux2_L[:, :, 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_L[:, :, n_nodes + 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_R[:, :, 1, :, element] .= zero(uEltype) - fluxes.antidiffusive_flux2_R[:, :, n_nodes + 1, :, element] .= zero(uEltype) - - fluxes.antidiffusive_flux3_L[:, :, :, 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux3_L[:, :, :, n_nodes + 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux3_R[:, :, :, 1, element] .= zero(uEltype) - fluxes.antidiffusive_flux3_R[:, :, :, n_nodes + 1, element] .= zero(uEltype) + return nothing +end + +function reset_antidiffusive_fluxes!(antidiffusive_flux1_L, antidiffusive_flux1_R, + antidiffusive_flux2_L, antidiffusive_flux2_R, + antidiffusive_flux3_L, antidiffusive_flux3_R) + uEltype = eltype(antidiffusive_flux1_L) + @threaded for element in axes(antidiffusive_flux1_L, 5) + antidiffusive_flux1_L[:, 1, :, :, element] .= zero(uEltype) + antidiffusive_flux1_L[:, end, :, :, element] .= zero(uEltype) + antidiffusive_flux1_R[:, 1, :, :, element] .= zero(uEltype) + antidiffusive_flux1_R[:, end, :, :, element] .= zero(uEltype) + + antidiffusive_flux2_L[:, :, 1, :, element] .= zero(uEltype) + antidiffusive_flux2_L[:, :, end, :, element] .= zero(uEltype) + antidiffusive_flux2_R[:, :, 1, :, element] .= zero(uEltype) + antidiffusive_flux2_R[:, :, end, :, element] .= zero(uEltype) + + antidiffusive_flux3_L[:, :, :, 1, element] .= zero(uEltype) + antidiffusive_flux3_L[:, :, :, end, element] .= zero(uEltype) + antidiffusive_flux3_R[:, :, :, 1, element] .= zero(uEltype) + antidiffusive_flux3_R[:, :, :, end, element] .= zero(uEltype) end return nothing diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index bccb96a1404..292a7952a40 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -24,7 +24,8 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, nnodes(dg), nnodes(dg)) for _ in 1:Threads.maxthreadid()] - antidiffusive_fluxes = ContainerAntidiffusiveFlux2D{uEltype}(0, + n_elements = nelements(cache_containers.elements) + antidiffusive_fluxes = ContainerAntidiffusiveFlux2D{uEltype}(n_elements, nvariables(equations), nnodes(dg)) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 31577758c26..9517ee04859 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -229,6 +229,8 @@ end function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{NDIMS}, basis::LobattoLegendreBasis, bound_keys) where {NDIMS} + # The number of elements is not yet known here. So, we initialize the container with 0 elements + # and resize it later while initializing the time integration method in `methods_SSP.jl`. subcell_limiter_coefficients = Trixi.ContainerSubcellLimiterIDP{NDIMS, real(basis)}(0, nnodes(basis), bound_keys) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 832140e4ff4..10df952532b 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -429,13 +429,13 @@ end ############################################################################### # Newton-bisection method -# 2D version @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, final_check, inverse_jacobian, dt, - equations, dg, cache, limiter) - (; inverse_weights) = dg.basis + equations::AbstractEquations{2}, + dg, cache, limiter) + (; inverse_weights) = dg.basis # Plays role of inverse DG-subcell sizes (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; gamma_constant_newton) = limiter diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index db7cf9021f8..ebe0e97ec34 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -132,9 +132,13 @@ function init(ode::ODEProblem, alg::SimpleAlgorithmSSP; kwargs...), false, true, false) - # resize container - resize!(integrator.p, integrator.p.solver.volume_integral, - nelements(integrator.p.solver, integrator.p.cache)) + # Resize container of volume integral for subcell limiting + _, _, dg, cache = mesh_equations_solver_cache(integrator.p) + if dg.volume_integral isa VolumeIntegralSubcellLimiting + # `subcell_limiter_coefficients` was created with 0 elements + resize!(dg.volume_integral.limiter.cache.subcell_limiter_coefficients, + nelements(dg, cache)) + end # Standard callbacks initialize_callbacks!(callback, integrator) @@ -258,11 +262,6 @@ function Base.resize!(integrator::SimpleIntegratorSSP, new_size) resize!(integrator.du, new_size) resize!(integrator.u_tmp, new_size) - # Resize container - # new_size = n_variables * n_nodes^n_dims * n_elements - n_elements = nelements(integrator.p.solver, integrator.p.cache) - resize!(integrator.p, integrator.p.solver.volume_integral, n_elements) - return nothing end end # @muladd