Closed bennibolm closed 1 month ago
This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.
NEWS.md
with its PR number.Created with :heart: by the Trixi.jl community.
Attention: Patch coverage is 99.05660%
with 2 lines
in your changes are missing coverage. Please review.
Project coverage is 96.11%. Comparing base (
c2513e2
) to head (83a8cef
).
Files | Patch % | Lines |
---|---|---|
...rc/solvers/dgsem_structured/subcell_limiters_2d.jl | 98.41% | 2 Missing :warning: |
:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.
Benchmarking if using dispatched routines is faster than calling inverse_jacobian
inside for node for loop.
Running kelving helmholtz on TreeMesh with local limiting
using Trixi
equations = CompressibleEulerEquations2D(1.4)
basis = LobattoLegendreBasis(3)
limiter_idp = SubcellLimiterIDP(equations, basis;
local_twosided_variables_cons = ["rho"],
positivity_variables_cons = ["rho"],
positivity_variables_nonlinear = [pressure])
trixi_include("../examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl", tspan = (0.0, 0.2), limiter_idp=limiter_idp)
Benchmarking:
limiter = semi.solver.volume_integral.limiter
alpha = limiter.cache.subcell_limiter_coefficients.alpha
u = Trixi.wrap_array(sol.u[end], mesh, equations, solver, semi.cache)
t=0.0
dt=0.1
variable = 1
using BenchmarkTools
@btime Trixi.idp_local_twosided!($alpha, $limiter, $u, $t, $dt, $semi, $mesh, $variable)
First version with
@threaded for element in eachelement(dg, semi.cache)
for j in eachnode(dg), i in eachnode(dg)
inverse_jacobian = cache.elements.inverse_jacobian[element]
idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache,
variable, var_min, var_max, i, j, element)
end
end
and
@threaded for element in eachelement(dg, semi.cache)
inverse_jacobian = cache.elements.inverse_jacobian[element]
for j in eachnode(dg), i in eachnode(dg)
idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache,
variable, var_min, var_max, i, j, element)
end
end
results in equivalent numbers: 338.607 μs (4 allocations: 136 bytes)
(The allocations are coming from creating the symbol in
var_min = variable_bounds[Symbol(variable_string, "_min")]
var_max = variable_bounds[Symbol(variable_string, "_max")]
Conclusion: Dispatching is not faster for the TreeMesh implementation. So, it's possible to use something like get_inverse_jacobian
to get the jacobian for all mesh types within the inner loop.
This saves a lot of lines in this PR.
Old description:
Getting the inverse Jacobian makes this PR a bit longer. That's because for TreeMesh, inverse_jacobian
is fixed within one element, while it is different for each node for StructuredMesh.
To avoid something like
if mesh isa TreeMesh
inverse_jacobian = cache.elements.inverse_jacobian[element]
end
for j in eachnode(dg), i in eachnode(dg)
if mesh isa StructuredMesh
inverse_jacobian = cache.elements.inverse_jacobian[i, j, element]
end
[...]
end
and also avoid unnecessary calls of inverse jacobian in TreeMesh simulations, I dispatched for the mesh type and extracted the code [...]
to an inner function (This is also done for e.g. the routine apply_jacobian!
).
UPDATE:
Due to https://github.com/trixi-framework/Trixi.jl/pull/1946#issuecomment-2112619960 and the addition of get_inverse_jacobian
, no dispatching is needed anymore.
This PR adds subcell limiting support for
Structuredmesh
. This includes the curved mesh version ofcalcflux_fhat
and the computation of local bounds at interfaces and boundaries (calc_bounds_twosided_interface
andcalc_bounds_onesided_interface
). Moreover, to get the inverse Jacobian for each node correctly, we useget_inverse_jacobian()
for all mesh types.The rest are just resorting and reformatting ugly lines of code.
For now, there is no support for non-conservative systems.