From 349373a61e84b62377967f69d23825973670683e Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 22 May 2025 12:34:35 +0100 Subject: [PATCH 01/38] Added steps to compute the gradient for hyperbolic systems. --- .../compressible_euler_multicomponent_2d.jl | 29 +++++++++++++++++++ .../compressible_navier_stokes_2d.jl | 4 +-- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index fdb42776a61..d8d80e0a8be 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -865,4 +865,33 @@ end v = u[orientation] / rho return v end + +@inline function enstrophy(u, equations::CompressibleEulerMulticomponentEquations2D) + # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v + + # Since gradients are calculated for parabolic equations we need to use + # some of their infrastructure. + viscous_container = init_viscous_container_2d(nvariables(equations), + nnodes(elements), nelements(elements), + uEltype) + @unpack u_transformed, gradients, flux_viscous = viscous_container + + calc_gradient!(gradients, u_transformed, 0.0, + mesh::TreeMesh{2}, equations, + boundary_conditions, dg::DG, parabolic_scheme, + cache, cache) + + omega = vorticity(u, gradients, equations) + + return 0.5f0 * u[1] * omega^2 +end + +@inline function vorticity(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) + # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) + + return dv2dx - dv1dy +end + end # @muladd diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index c5ca6586bfa..5f595dfdbe6 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -288,14 +288,14 @@ end return SVector(v1, v2) end -@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D) +@inline function enstrophy(u, equations::CompressibleNavierStokesDiffusion2D; gradients) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations) return 0.5f0 * u[1] * omega^2 end -@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D) +@inline function vorticity(u, equations::CompressibleNavierStokesDiffusion2D; gradients) # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) From 9277c5b76822aca55656dca9a0e254e95b06739a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 22 May 2025 13:49:44 +0100 Subject: [PATCH 02/38] Moved potential gradient calculations into integration routine. --- src/callbacks_step/analysis_dg2d.jl | 28 +++++++++++++++---- .../compressible_euler_multicomponent_2d.jl | 14 +--------- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index dd166bbe389..6c47f338dd4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -246,11 +246,29 @@ function integrate(func::Func, u, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} - integrate_via_indices(u, mesh, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg - u_local = get_node_vars(u, equations, dg, i, j, element) - return func(u_local, equations) - end + @autoinfiltrate + if length(m.sig.parameters) == 2 + integrate_via_indices(u, mesh, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, element) + return func(u_local, equations) + end + if length(m.sig.parameters) == 3 + integrate_via_indices(u, mesh, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, element) + # Since gradients are calculated for parabolic equations we need to use + # some of their infrastructure. + viscous_container = init_viscous_container_2d(nvariables(equations), + nnodes(cache.elements), nelements(cache.elements), + eltype(cache.elements)) + @unpack u_transformed, gradients, flux_viscous = viscous_container + calc_gradient!(gradients, u_transformed, 0.0, + mesh::TreeMesh{2}, equations, + boundary_conditions, dg::DG, parabolic_scheme, + cache, cache) + return func(u_local, gradients, equations) + end end function integrate(func::Func, u, diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index d8d80e0a8be..1d54705cdac 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -866,21 +866,9 @@ end return v end -@inline function enstrophy(u, equations::CompressibleEulerMulticomponentEquations2D) +@inline function enstrophy(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v - # Since gradients are calculated for parabolic equations we need to use - # some of their infrastructure. - viscous_container = init_viscous_container_2d(nvariables(equations), - nnodes(elements), nelements(elements), - uEltype) - @unpack u_transformed, gradients, flux_viscous = viscous_container - - calc_gradient!(gradients, u_transformed, 0.0, - mesh::TreeMesh{2}, equations, - boundary_conditions, dg::DG, parabolic_scheme, - cache, cache) - omega = vorticity(u, gradients, equations) return 0.5f0 * u[1] * omega^2 From 9826f54401b53b3cc6cf55ff46fb0da4ab98789b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 27 May 2025 15:51:12 +0100 Subject: [PATCH 03/38] Distinguish between call using functions with 2 or 3 parameters. --- src/callbacks_step/analysis_dg2d.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 6c47f338dd4..ab48d51a507 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -247,13 +247,15 @@ function integrate(func::Func, u, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} @autoinfiltrate - if length(m.sig.parameters) == 2 + m = methods(func) + if length(m[1].sig.parameters) == 2 integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) end - if length(m.sig.parameters) == 3 + end + if length(m[1].sig.parameters) == 3 integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) @@ -264,11 +266,12 @@ function integrate(func::Func, u, eltype(cache.elements)) @unpack u_transformed, gradients, flux_viscous = viscous_container calc_gradient!(gradients, u_transformed, 0.0, - mesh::TreeMesh{2}, equations, - boundary_conditions, dg::DG, parabolic_scheme, - cache, cache) + mesh::TreeMesh{2}, equations, + boundary_conditions, dg::DG, parabolic_scheme, + cache, cache) return func(u_local, gradients, equations) end + end end function integrate(func::Func, u, From 5c43164a387b6a6d8d9759fc4cd1479b433b4ae3 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Sat, 31 May 2025 00:25:17 +0100 Subject: [PATCH 04/38] Reverted function signature. --- src/equations/compressible_navier_stokes_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 5f595dfdbe6..c5ca6586bfa 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -288,14 +288,14 @@ end return SVector(v1, v2) end -@inline function enstrophy(u, equations::CompressibleNavierStokesDiffusion2D; gradients) +@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations) return 0.5f0 * u[1] * omega^2 end -@inline function vorticity(u, equations::CompressibleNavierStokesDiffusion2D; gradients) +@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D) # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) From c889484c56bf3ad36a19963e43657a4d81600e19 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 2 Jun 2025 14:57:01 +0100 Subject: [PATCH 05/38] Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/compressible_euler_multicomponent_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 1d54705cdac..e1d8726e081 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -866,7 +866,8 @@ end return v end -@inline function enstrophy(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) +@inline function enstrophy(u, gradients, + equations::CompressibleEulerMulticomponentEquations2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations) From 2821bf026ab16f590f03f31ee345f3e6957dec43 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 2 Jun 2025 14:57:07 +0100 Subject: [PATCH 06/38] Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/compressible_euler_multicomponent_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index e1d8726e081..3134114ad4f 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -875,7 +875,8 @@ end return 0.5f0 * u[1] * omega^2 end -@inline function vorticity(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) +@inline function vorticity(u, gradients, + equations::CompressibleEulerMulticomponentEquations2D) # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) From 33ead0dc893a57bb64cb1ace399c05dd69546906 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 2 Jun 2025 14:57:13 +0100 Subject: [PATCH 07/38] Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/compressible_euler_multicomponent_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 3134114ad4f..dceb8f429fd 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -883,5 +883,4 @@ end return dv2dx - dv1dy end - end # @muladd From 818c698ffb244529d257c89c9c77960c9285df7e Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 2 Jun 2025 16:09:32 +0100 Subject: [PATCH 08/38] Added Andrew's Fortran code for derivatives. Attempt of using existing infrastructure for parabolic terms to compute derivatives. --- src/callbacks_step/analysis_dg2d.jl | 41 +++++++++++++++++++- src/semidiscretization/semidiscretization.jl | 2 +- 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index ab48d51a507..6b06cd2a0d6 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -245,8 +245,10 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG, cache; normalize = true) where {Func} + equations, dg::DG, cache; normalize = true, semi) where {Func} @autoinfiltrate + @unpack solver = semi + @unpack boundaries = cache m = methods(func) if length(m[1].sig.parameters) == 2 integrate_via_indices(u, mesh, equations, dg, cache; @@ -267,13 +269,48 @@ function integrate(func::Func, u, @unpack u_transformed, gradients, flux_viscous = viscous_container calc_gradient!(gradients, u_transformed, 0.0, mesh::TreeMesh{2}, equations, - boundary_conditions, dg::DG, parabolic_scheme, + boundaries, dg::DG, solver, cache, cache) return func(u_local, gradients, equations) end end end +# Andrew's functions for computing the derivatives. +# function DGSpaceDerivative_WeakForm!( +# F_prime::Matrix{Float64}, # size: (N+1, nEqn) +# FL::Vector{Float64}, # size: nEqn +# FR::Vector{Float64}, # size: nEqn +# Flux::Matrix{Float64}, # size: (N+1, nEqn) +# D::Matrix{Float64}, # size: (N+1, N+1) +# l_minus::Vector{Float64}, # size: N+1 +# l_plus::Vector{Float64}, # size: N+1 +# weights::Vector{Float64}, # size: N+1 +# nEqn::Int, +# N::Int +# ) +# # Volume derivative term using D matrix +# for k in 1:nEqn +# MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) +# end +# +# # Surface contribution from numerical fluxes +# for j in 0:N +# for k in 1:nEqn +# F_prime[j+1, k] += (FR[k] * l_plus[j+1] + FL[k] * l_minus[j+1]) / weights[j+1] +# end +# end +# end +# +# function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Matrix{Float64}, N::Int) +# for i in 0:N +# t = 0.0 +# for j in 0:N +# t += D[i+1, j+1] * Phi[j+1] # Adjust for 1-based indexing +# end +# Phi_prime[i+1] = t +# end + function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, equations_parabolic, diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index f41c7ea4a7f..75174d6c076 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -58,7 +58,7 @@ function integrate(func::Func, u_ode, semi::AbstractSemidiscretization; mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - integrate(func, u, mesh, equations, solver, cache, normalize = normalize) + integrate(func, u, mesh, equations, solver, cache, normalize = normalize; semi=semi) end function integrate(u_ode, semi::AbstractSemidiscretization; normalize = true) From 103c918423092175ee228d0b2324eabc5ea7e4f2 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 3 Jun 2025 15:53:25 +0100 Subject: [PATCH 09/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 6b06cd2a0d6..71cceb6c7c4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -252,7 +252,7 @@ function integrate(func::Func, u, m = methods(func) if length(m[1].sig.parameters) == 2 integrate_via_indices(u, mesh, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg + normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) end From 76adb1050820e7f28262e54c1c37d6d17acdc2f0 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 3 Jun 2025 15:53:41 +0100 Subject: [PATCH 10/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 71cceb6c7c4..a15ed3c1073 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -259,7 +259,7 @@ function integrate(func::Func, u, end if length(m[1].sig.parameters) == 3 integrate_via_indices(u, mesh, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg + normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) # Since gradients are calculated for parabolic equations we need to use # some of their infrastructure. From 373c1d64a29e480490b4b0f120c8c1bdfa43e285 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 3 Jun 2025 15:53:52 +0100 Subject: [PATCH 11/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index a15ed3c1073..12a8898eec4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -264,8 +264,9 @@ function integrate(func::Func, u, # Since gradients are calculated for parabolic equations we need to use # some of their infrastructure. viscous_container = init_viscous_container_2d(nvariables(equations), - nnodes(cache.elements), nelements(cache.elements), - eltype(cache.elements)) + nnodes(cache.elements), + nelements(cache.elements), + eltype(cache.elements)) @unpack u_transformed, gradients, flux_viscous = viscous_container calc_gradient!(gradients, u_transformed, 0.0, mesh::TreeMesh{2}, equations, From 4e6245b304e084ac39f2ce6c5dd19866ee8a7c24 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 3 Jun 2025 15:57:06 +0100 Subject: [PATCH 12/38] Update src/semidiscretization/semidiscretization.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/semidiscretization/semidiscretization.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 75174d6c076..d575aa6d1da 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -58,7 +58,8 @@ function integrate(func::Func, u_ode, semi::AbstractSemidiscretization; mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - integrate(func, u, mesh, equations, solver, cache, normalize = normalize; semi=semi) + integrate(func, u, mesh, equations, solver, cache, normalize = normalize; + semi = semi) end function integrate(u_ode, semi::AbstractSemidiscretization; normalize = true) From a93e340030803d21cfeb89df503696a7da2b3d5d Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 4 Jun 2025 17:19:38 +0100 Subject: [PATCH 13/38] Sterted conversion of Andrew's derivatives method. --- src/callbacks_step/analysis_dg2d.jl | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 12a8898eec4..1635b24341f 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -277,19 +277,25 @@ function integrate(func::Func, u, end end -# Andrew's functions for computing the derivatives. +# # Andrew's functions for computing the derivatives. # function DGSpaceDerivative_WeakForm!( -# F_prime::Matrix{Float64}, # size: (N+1, nEqn) +# # F_prime::Matrix{Float64}, # size: (N+1, nEqn) # FL::Vector{Float64}, # size: nEqn # FR::Vector{Float64}, # size: nEqn -# Flux::Matrix{Float64}, # size: (N+1, nEqn) -# D::Matrix{Float64}, # size: (N+1, N+1) +# # Flux::Matrix{Float64}, # size: (N+1, nEqn) +# # D::Matrix{Float64}, # size: (N+1, N+1) # l_minus::Vector{Float64}, # size: N+1 # l_plus::Vector{Float64}, # size: N+1 # weights::Vector{Float64}, # size: N+1 # nEqn::Int, # N::Int # ) +# # Translations. +# D = basis.derivative_dhat +# J = elements.jacobian_matrix +# Flux = +# F_prime = D*Flux +# # # Volume derivative term using D matrix # for k in 1:nEqn # MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) @@ -311,6 +317,7 @@ end # end # Phi_prime[i+1] = t # end +# end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, From d7f3f55426a9053c662258e7010b7ec1154fb6ea Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 11 Jun 2025 10:19:07 +0100 Subject: [PATCH 14/38] Fuirther implementations of Andrew's derivative funciton. --- src/callbacks_step/analysis_dg2d.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 1635b24341f..986b7167d60 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -279,22 +279,23 @@ end # # Andrew's functions for computing the derivatives. # function DGSpaceDerivative_WeakForm!( -# # F_prime::Matrix{Float64}, # size: (N+1, nEqn) # FL::Vector{Float64}, # size: nEqn # FR::Vector{Float64}, # size: nEqn -# # Flux::Matrix{Float64}, # size: (N+1, nEqn) -# # D::Matrix{Float64}, # size: (N+1, N+1) # l_minus::Vector{Float64}, # size: N+1 # l_plus::Vector{Float64}, # size: N+1 -# weights::Vector{Float64}, # size: N+1 # nEqn::Int, # N::Int # ) +# # Get the required variables. +# @unpack derivative_dhat = dg.basis +# # # Translations. -# D = basis.derivative_dhat +# D = derivative_dhat # J = elements.jacobian_matrix # Flux = # F_prime = D*Flux +# weights = dg.basis.weights +# spA_Dhat = dg.basis.derivative_dhat # # # Volume derivative term using D matrix # for k in 1:nEqn From 5c4e4db36aa0e6c5fe74a4d7d5dd1b59e03887e4 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 12 Jun 2025 13:19:44 +0100 Subject: [PATCH 15/38] Further translations of Fortran variables to Trixi.jl variables. --- src/callbacks_step/analysis_dg2d.jl | 88 +++++++++++++++-------------- 1 file changed, 46 insertions(+), 42 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 986b7167d60..d6c76f53b3b 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -277,48 +277,52 @@ function integrate(func::Func, u, end end -# # Andrew's functions for computing the derivatives. -# function DGSpaceDerivative_WeakForm!( -# FL::Vector{Float64}, # size: nEqn -# FR::Vector{Float64}, # size: nEqn -# l_minus::Vector{Float64}, # size: N+1 -# l_plus::Vector{Float64}, # size: N+1 -# nEqn::Int, -# N::Int -# ) -# # Get the required variables. -# @unpack derivative_dhat = dg.basis -# -# # Translations. -# D = derivative_dhat -# J = elements.jacobian_matrix -# Flux = -# F_prime = D*Flux -# weights = dg.basis.weights -# spA_Dhat = dg.basis.derivative_dhat -# -# # Volume derivative term using D matrix -# for k in 1:nEqn -# MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) -# end -# -# # Surface contribution from numerical fluxes -# for j in 0:N -# for k in 1:nEqn -# F_prime[j+1, k] += (FR[k] * l_plus[j+1] + FL[k] * l_minus[j+1]) / weights[j+1] -# end -# end -# end -# -# function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Matrix{Float64}, N::Int) -# for i in 0:N -# t = 0.0 -# for j in 0:N -# t += D[i+1, j+1] * Phi[j+1] # Adjust for 1-based indexing -# end -# Phi_prime[i+1] = t -# end -# end +# Andrew's functions for computing the derivatives. +function DGSpaceDerivative_WeakForm!( + FL::Vector{Float64}, # size: nEqn + FR::Vector{Float64}, # size: nEqn + l_minus::Vector{Float64}, # size: N+1 + l_plus::Vector{Float64}, # size: N+1 + nEqn::Int, + N::Int +) + # Get the required variables. + @unpack derivative_dhat = dg.basis + + # Translations. + D = derivative_dhat + spA_Dhat = D + J = elements.jacobian_matrix + Flux = dg.volume_integral.volume_flux_dg(???) + F_prime = spA_Dhat * Flux + weights = dg.basis.weights + geom%X_xi = elements.contravariant_vectors ??? + l_minus = dg.basis.basis_at_ξ̂_min + l_plus = dg.basis.basis_at_ξ̂_max + Fstarb = + + # Volume derivative term using D matrix + for k in 1:nEqn + MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) + end + + # Surface contribution from numerical fluxes + for j in 0:N + for k in 1:nEqn + F_prime[j+1, k] += (FR[k] * l_plus[j+1] + FL[k] * l_minus[j+1]) / weights[j+1] + end + end +end + +function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Matrix{Float64}, N::Int) + for i in 0:N + t = 0.0 + for j in 0:N + t += D[i+1, j+1] * Phi[j+1] # Adjust for 1-based indexing + end + Phi_prime[i+1] = t + end +end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, From 00d61c9fba6d8cb8866457f8b3ba99ed80caad03 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Jun 2025 15:47:03 +0100 Subject: [PATCH 16/38] Added calculation of fluxes by direciton. --- src/callbacks_step/analysis_dg2d.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index d6c76f53b3b..f92a500530e 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -285,6 +285,7 @@ function DGSpaceDerivative_WeakForm!( l_plus::Vector{Float64}, # size: N+1 nEqn::Int, N::Int + direction::Int ) # Get the required variables. @unpack derivative_dhat = dg.basis @@ -292,14 +293,12 @@ function DGSpaceDerivative_WeakForm!( # Translations. D = derivative_dhat spA_Dhat = D - J = elements.jacobian_matrix - Flux = dg.volume_integral.volume_flux_dg(???) + J = 1./cache.elements.inverse_jacobian +# Flux = dg.volume_integral.volume_flux_dg(???) F_prime = spA_Dhat * Flux weights = dg.basis.weights - geom%X_xi = elements.contravariant_vectors ??? - l_minus = dg.basis.basis_at_ξ̂_min - l_plus = dg.basis.basis_at_ξ̂_max - Fstarb = + l_minus = l_minus = dg.basis.boundary_interpolation[:, 1] + l_plus = dg.basis.boundary_interpolation[:, 2] # Volume derivative term using D matrix for k in 1:nEqn @@ -324,6 +323,13 @@ function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Mat end end +function compute_flux_array!(Flux, u, direction, equations) + Np, nvars = size(Flux) + for i in 1:Np + Flux[i, :] .= flux(u[i, :], 1, direction, equations) + end +end + function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, equations_parabolic, From 142493058b87471877017eca62978fa7f478d25c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 17 Jun 2025 14:39:07 +0100 Subject: [PATCH 17/38] Corrected calculation of gradients. --- src/callbacks_step/analysis_dg2d.jl | 82 ++++++++++++++++++----------- 1 file changed, 52 insertions(+), 30 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index f92a500530e..f57f4a62ec1 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -246,7 +246,6 @@ function integrate(func::Func, u, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true, semi) where {Func} - @autoinfiltrate @unpack solver = semi @unpack boundaries = cache m = methods(func) @@ -261,17 +260,9 @@ function integrate(func::Func, u, integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) - # Since gradients are calculated for parabolic equations we need to use - # some of their infrastructure. - viscous_container = init_viscous_container_2d(nvariables(equations), - nnodes(cache.elements), - nelements(cache.elements), - eltype(cache.elements)) - @unpack u_transformed, gradients, flux_viscous = viscous_container - calc_gradient!(gradients, u_transformed, 0.0, - mesh::TreeMesh{2}, equations, - boundaries, dg::DG, solver, - cache, cache) + @autoinfiltrate + gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) + return func(u_local, gradients, equations) end end @@ -279,13 +270,11 @@ end # Andrew's functions for computing the derivatives. function DGSpaceDerivative_WeakForm!( - FL::Vector{Float64}, # size: nEqn - FR::Vector{Float64}, # size: nEqn - l_minus::Vector{Float64}, # size: N+1 - l_plus::Vector{Float64}, # size: N+1 - nEqn::Int, - N::Int - direction::Int + dg, + cache, + u, + direction::Int, + equations ) # Get the required variables. @unpack derivative_dhat = dg.basis @@ -293,24 +282,57 @@ function DGSpaceDerivative_WeakForm!( # Translations. D = derivative_dhat spA_Dhat = D - J = 1./cache.elements.inverse_jacobian -# Flux = dg.volume_integral.volume_flux_dg(???) - F_prime = spA_Dhat * Flux + J = 1 ./ cache.elements.inverse_jacobian weights = dg.basis.weights - l_minus = l_minus = dg.basis.boundary_interpolation[:, 1] + l_minus = dg.basis.boundary_interpolation[:, 1] l_plus = dg.basis.boundary_interpolation[:, 2] - # Volume derivative term using D matrix - for k in 1:nEqn - MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) + Np, nvars = size(u) + F = similar(u) + F_prime = similar(u) + + # Multiply Dhat * F for each variable + for k in 1:nvars + for i in 1:Np + F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) + end end - # Surface contribution from numerical fluxes - for j in 0:N - for k in 1:nEqn - F_prime[j+1, k] += (FR[k] * l_plus[j+1] + FL[k] * l_minus[j+1]) / weights[j+1] + # Store into gradients array + for i in 1:Np + for k in 1:nvars + gradients[k, direction, i] = F_prime[i, k] end end + + + +# # Evaluate flux at each point +# for i in 1:Np +# F[i, :] .= flux_func(u[i, :], orientation, direction, equations) +# end +# +# Flux = compute_flux_array!(Flux, u, direction, equations) +# # Volume derivative term using D matrix +# for k in 1:nEqn +# MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) +# end +# +# # Multiply Dhat * F for each variable +# F_prime = spA_Dhat * Flux +# for k in 1:nvars +# for i in 1:Np +# F_prime[i, k] = sum(Dhat[i, j] * F[j, k] for j in 1:Np) +# end +# end +#= + for j in 1:Np + for k in 1:nvars + gradients[k, direction, j] = F_prime[j, k] + end + end=# + + return gradients end function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Matrix{Float64}, N::Int) From 0dcc5ee60f9cc81011b547babfddfb352e62feae Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 23 Jun 2025 16:37:36 +0100 Subject: [PATCH 18/38] Corrected the gradient function. --- src/callbacks_step/analysis_dg2d.jl | 60 +++++++++++------------------ 1 file changed, 22 insertions(+), 38 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index f57f4a62ec1..2501841c9ce 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -260,9 +260,9 @@ function integrate(func::Func, u, integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) - @autoinfiltrate gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) + @autoinfiltrate return func(u_local, gradients, equations) end end @@ -287,50 +287,34 @@ function DGSpaceDerivative_WeakForm!( l_minus = dg.basis.boundary_interpolation[:, 1] l_plus = dg.basis.boundary_interpolation[:, 2] - Np, nvars = size(u) - F = similar(u) - F_prime = similar(u) - - # Multiply Dhat * F for each variable - for k in 1:nvars - for i in 1:Np - F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) - end - end - - # Store into gradients array - for i in 1:Np - for k in 1:nvars - gradients[k, direction, i] = F_prime[i, k] - end - end - - - -# # Evaluate flux at each point -# for i in 1:Np -# F[i, :] .= flux_func(u[i, :], orientation, direction, equations) -# end -# -# Flux = compute_flux_array!(Flux, u, direction, equations) -# # Volume derivative term using D matrix -# for k in 1:nEqn -# MxVDerivative!(F_prime[:, k], Flux[:, k], D, N) -# end +# nvars, Np = size(u) +# @autoinfiltrate +# F = similar(u) +# F_prime = similar(u) # # # Multiply Dhat * F for each variable -# F_prime = spA_Dhat * Flux # for k in 1:nvars # for i in 1:Np -# F_prime[i, k] = sum(Dhat[i, j] * F[j, k] for j in 1:Np) +# F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) # end # end -#= - for j in 1:Np - for k in 1:nvars - gradients[k, direction, j] = F_prime[j, k] +# +# # Store into gradients array +# for i in 1:Np +# for k in 1:nvars +# gradients[k, i] = F_prime[i, k] +# end +# end + + n_vars, Np, _, n_elements = size(u) + gradients = similar(u) + + for elem in 1:n_elements + for v in 1:n_vars + # Contract D in the ξ̂ (first spatial) direction + gradients[v, :, :, elem] .= D * u[v, :, :, elem] end - end=# + end return gradients end From 470c96db0a8d0ae201f097c79c42c519d5902326 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 26 Jun 2025 11:29:33 +0100 Subject: [PATCH 19/38] Added 'semi' option to function call. --- src/callbacks_step/analysis.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index bbc2ac4c8fe..7d9e13d1435 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -617,10 +617,10 @@ end # Trixi.analyze, Trixi.pretty_form_utf, Trixi.pretty_form_ascii function analyze(quantity, du, u, t, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - analyze(quantity, du, u, t, mesh, equations, solver, cache) + analyze(quantity, du, u, t, mesh, equations, solver, cache; semi=semi) end -function analyze(quantity, du, u, t, mesh, equations, solver, cache) - integrate(quantity, u, mesh, equations, solver, cache, normalize = true) +function analyze(quantity, du, u, t, mesh, equations, solver, cache; semi=nothing) + integrate(quantity, u, mesh, equations, solver, cache, normalize = true, semi=semi) end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) From db9b07d01ce083c0138cd7c13d555ef2b57f24e7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 26 Jun 2025 11:29:53 +0100 Subject: [PATCH 20/38] Take care of cases when function is cons2cons. --- src/callbacks_step/analysis_dg2d.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 2501841c9ce..f07bf0be661 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -249,21 +249,22 @@ function integrate(func::Func, u, @unpack solver = semi @unpack boundaries = cache m = methods(func) - if length(m[1].sig.parameters) == 2 - integrate_via_indices(u, mesh, equations, dg, cache; + @autoinfiltrate + if (m[1].nargs == 2) || (func == cons2cons) + return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) - return func(u_local, equations) + + func(u_local, equations) end end - if length(m[1].sig.parameters) == 3 - integrate_via_indices(u, mesh, equations, dg, cache; + if (m[1].nargs == 3) && (func != cons2cons) + return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) - @autoinfiltrate - return func(u_local, gradients, equations) + func(u_local, gradients, equations) end end end From a2e26c29a8cac9ab473d5859ab941aca097baaa7 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 30 Jun 2025 14:19:06 +0100 Subject: [PATCH 21/38] Update src/callbacks_step/analysis.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 7d9e13d1435..42c36d06dfb 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -617,7 +617,7 @@ end # Trixi.analyze, Trixi.pretty_form_utf, Trixi.pretty_form_ascii function analyze(quantity, du, u, t, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - analyze(quantity, du, u, t, mesh, equations, solver, cache; semi=semi) + analyze(quantity, du, u, t, mesh, equations, solver, cache; semi = semi) end function analyze(quantity, du, u, t, mesh, equations, solver, cache; semi=nothing) integrate(quantity, u, mesh, equations, solver, cache, normalize = true, semi=semi) From 71f7e7af990eee849c4e777b3e143f9f2434a21f Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 30 Jun 2025 14:19:14 +0100 Subject: [PATCH 22/38] Update src/callbacks_step/analysis.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 42c36d06dfb..9550428ed05 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -619,8 +619,9 @@ function analyze(quantity, du, u, t, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) analyze(quantity, du, u, t, mesh, equations, solver, cache; semi = semi) end -function analyze(quantity, du, u, t, mesh, equations, solver, cache; semi=nothing) - integrate(quantity, u, mesh, equations, solver, cache, normalize = true, semi=semi) +function analyze(quantity, du, u, t, mesh, equations, solver, cache; semi = nothing) + integrate(quantity, u, mesh, equations, solver, cache, normalize = true, + semi = semi) end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) From c705cdf1c164888de7212f6dc0de74a48ec9fd42 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 30 Jun 2025 14:19:23 +0100 Subject: [PATCH 23/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index f07bf0be661..b5adb538d29 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -252,7 +252,8 @@ function integrate(func::Func, u, @autoinfiltrate if (m[1].nargs == 2) || (func == cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg + normalize = normalize) do u, i, j, element, + equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) func(u_local, equations) From 0c50c08e54fdb0f0904af71d27f747b4aecb26bb Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Jul 2025 12:36:39 +0100 Subject: [PATCH 24/38] Removed redundant 'semi' parameter. --- src/callbacks_step/analysis.jl | 7 +++---- src/callbacks_step/analysis_dg2d.jl | 4 +--- src/semidiscretization/semidiscretization.jl | 3 +-- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 9550428ed05..bbc2ac4c8fe 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -617,11 +617,10 @@ end # Trixi.analyze, Trixi.pretty_form_utf, Trixi.pretty_form_ascii function analyze(quantity, du, u, t, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - analyze(quantity, du, u, t, mesh, equations, solver, cache; semi = semi) + analyze(quantity, du, u, t, mesh, equations, solver, cache) end -function analyze(quantity, du, u, t, mesh, equations, solver, cache; semi = nothing) - integrate(quantity, u, mesh, equations, solver, cache, normalize = true, - semi = semi) +function analyze(quantity, du, u, t, mesh, equations, solver, cache) + integrate(quantity, u, mesh, equations, solver, cache, normalize = true) end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index b5adb538d29..d338abeb938 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -245,11 +245,9 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG, cache; normalize = true, semi) where {Func} - @unpack solver = semi + equations, dg::DG, cache; normalize = true) where {Func} @unpack boundaries = cache m = methods(func) - @autoinfiltrate if (m[1].nargs == 2) || (func == cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index d575aa6d1da..f41c7ea4a7f 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -58,8 +58,7 @@ function integrate(func::Func, u_ode, semi::AbstractSemidiscretization; mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - integrate(func, u, mesh, equations, solver, cache, normalize = normalize; - semi = semi) + integrate(func, u, mesh, equations, solver, cache, normalize = normalize) end function integrate(u_ode, semi::AbstractSemidiscretization; normalize = true) From 96588746de6728dfa7b544bdbc775af1eb977353 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:45:13 +0100 Subject: [PATCH 25/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index d338abeb938..7f20fd8f3d5 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -259,7 +259,8 @@ function integrate(func::Func, u, end if (m[1].nargs == 3) && (func != cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg + normalize = normalize) do u, i, j, element, + equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) From 423dd30553398a0939e3ffc5384a0e968f787c43 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:45:20 +0100 Subject: [PATCH 26/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 7f20fd8f3d5..3471c98044f 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -270,13 +270,11 @@ function integrate(func::Func, u, end # Andrew's functions for computing the derivatives. -function DGSpaceDerivative_WeakForm!( - dg, - cache, - u, - direction::Int, - equations -) +function DGSpaceDerivative_WeakForm!(dg, + cache, + u, + direction::Int, + equations) # Get the required variables. @unpack derivative_dhat = dg.basis From 1081df072c15ad730202dbd60bf9316f2d0bc8af Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:45:43 +0100 Subject: [PATCH 27/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 36 ++++++++++++++--------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 3471c98044f..a123741586a 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -286,24 +286,24 @@ function DGSpaceDerivative_WeakForm!(dg, l_minus = dg.basis.boundary_interpolation[:, 1] l_plus = dg.basis.boundary_interpolation[:, 2] -# nvars, Np = size(u) -# @autoinfiltrate -# F = similar(u) -# F_prime = similar(u) -# -# # Multiply Dhat * F for each variable -# for k in 1:nvars -# for i in 1:Np -# F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) -# end -# end -# -# # Store into gradients array -# for i in 1:Np -# for k in 1:nvars -# gradients[k, i] = F_prime[i, k] -# end -# end + # nvars, Np = size(u) + # @autoinfiltrate + # F = similar(u) + # F_prime = similar(u) + # + # # Multiply Dhat * F for each variable + # for k in 1:nvars + # for i in 1:Np + # F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) + # end + # end + # + # # Store into gradients array + # for i in 1:Np + # for k in 1:nvars + # gradients[k, i] = F_prime[i, k] + # end + # end n_vars, Np, _, n_elements = size(u) gradients = similar(u) From 6e31ee33416c31bee60ead64a305d954a7e0a23b Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:45:51 +0100 Subject: [PATCH 28/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index a123741586a..cc66c069e60 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -318,7 +318,8 @@ function DGSpaceDerivative_WeakForm!(dg, return gradients end -function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, D::Matrix{Float64}, N::Int) +function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, + D::Matrix{Float64}, N::Int) for i in 0:N t = 0.0 for j in 0:N From 3612f4ba7471a5ebb697e028d65970e6ae13f8c3 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:46:02 +0100 Subject: [PATCH 29/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index cc66c069e60..f4d24af6995 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -323,7 +323,7 @@ function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, for i in 0:N t = 0.0 for j in 0:N - t += D[i+1, j+1] * Phi[j+1] # Adjust for 1-based indexing + t += D[i + 1, j + 1] * Phi[j + 1] # Adjust for 1-based indexing end Phi_prime[i+1] = t end From 71ca315f816a48d8f87b6fe92a112553893965b4 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 3 Jul 2025 15:46:10 +0100 Subject: [PATCH 30/38] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index f4d24af6995..bb6c8cb0bd0 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -325,7 +325,7 @@ function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, for j in 0:N t += D[i + 1, j + 1] * Phi[j + 1] # Adjust for 1-based indexing end - Phi_prime[i+1] = t + Phi_prime[i + 1] = t end end From b8a2716382982ab5d08db1a590c0173ca96c6813 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 3 Jul 2025 15:50:38 +0100 Subject: [PATCH 31/38] Removed redundant code. --- src/callbacks_step/analysis_dg2d.jl | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index d338abeb938..48c1a677bf9 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -287,25 +287,6 @@ function DGSpaceDerivative_WeakForm!( l_minus = dg.basis.boundary_interpolation[:, 1] l_plus = dg.basis.boundary_interpolation[:, 2] -# nvars, Np = size(u) -# @autoinfiltrate -# F = similar(u) -# F_prime = similar(u) -# -# # Multiply Dhat * F for each variable -# for k in 1:nvars -# for i in 1:Np -# F_prime[i, k] = sum(D[i, j] * F[j, k] for j in 1:Np) -# end -# end -# -# # Store into gradients array -# for i in 1:Np -# for k in 1:nvars -# gradients[k, i] = F_prime[i, k] -# end -# end - n_vars, Np, _, n_elements = size(u) gradients = similar(u) From 2d65ccc6598f67737e3e4308519de0e74234aa5f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 4 Jul 2025 08:57:07 +0100 Subject: [PATCH 32/38] Remoevd some redundant code. --- src/Trixi.jl | 2 +- src/callbacks_step/analysis_dg2d.jl | 27 ++++++++++--------- .../compressible_euler_multicomponent_2d.jl | 2 +- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 8f13835dbae..e5a46172f49 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -235,7 +235,7 @@ export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight, waterheight_pressure export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, - enstrophy, magnetic_field, divergence_cleaning_field + enstrophy, magnetic_field, divergence_cleaning_field, enstrophy_multi_euler export lake_at_rest_error export ncomponents, eachcomponent diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 038fbb85e67..9afc978b4e4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -247,6 +247,7 @@ function integrate(func::Func, u, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} @unpack boundaries = cache + @autoinfiltrate m = methods(func) if (m[1].nargs == 2) || (func == cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; @@ -257,14 +258,15 @@ function integrate(func::Func, u, func(u_local, equations) end end - if (m[1].nargs == 3) && (func != cons2cons) + if (m[1].nargs == 4) && (func != cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) - func(u_local, gradients, equations) +# func(u_local, gradients, equations) + return 0.0 end end end @@ -289,6 +291,7 @@ function DGSpaceDerivative_WeakForm!(dg, n_vars, Np, _, n_elements = size(u) gradients = similar(u) + @autoinfiltrate for elem in 1:n_elements for v in 1:n_vars # Contract D in the ξ̂ (first spatial) direction @@ -299,16 +302,16 @@ function DGSpaceDerivative_WeakForm!(dg, return gradients end -function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, - D::Matrix{Float64}, N::Int) - for i in 0:N - t = 0.0 - for j in 0:N - t += D[i + 1, j + 1] * Phi[j + 1] # Adjust for 1-based indexing - end - Phi_prime[i + 1] = t - end -end +# function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, +# D::Matrix{Float64}, N::Int) +# for i in 0:N +# t = 0.0 +# for j in 0:N +# t += D[i + 1, j + 1] * Phi[j + 1] # Adjust for 1-based indexing +# end +# Phi_prime[i + 1] = t +# end +# end function compute_flux_array!(Flux, u, direction, equations) Np, nvars = size(Flux) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index dceb8f429fd..844ea6e4e33 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -866,7 +866,7 @@ end return v end -@inline function enstrophy(u, gradients, +@inline function enstrophy_multi_euler(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v From 89cbfc06a18d332875c3584ef2882a4ccb85df09 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 10 Jul 2025 17:33:30 +0100 Subject: [PATCH 33/38] Compute now also the y-derivative. --- src/callbacks_step/analysis_dg2d.jl | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 9afc978b4e4..202ab5c7ab8 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -247,7 +247,6 @@ function integrate(func::Func, u, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} @unpack boundaries = cache - @autoinfiltrate m = methods(func) if (m[1].nargs == 2) || (func == cons2cons) return integrate_via_indices(u, mesh, equations, dg, cache; @@ -259,13 +258,16 @@ function integrate(func::Func, u, end end if (m[1].nargs == 4) && (func != cons2cons) + gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) + gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations) + gradients = Vector([gradients_x, gradients_y]) + @autoinfiltrate return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) - gradients = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) -# func(u_local, gradients, equations) + func(u_local, gradients, equations) return 0.0 end end @@ -291,7 +293,6 @@ function DGSpaceDerivative_WeakForm!(dg, n_vars, Np, _, n_elements = size(u) gradients = similar(u) - @autoinfiltrate for elem in 1:n_elements for v in 1:n_vars # Contract D in the ξ̂ (first spatial) direction @@ -302,17 +303,6 @@ function DGSpaceDerivative_WeakForm!(dg, return gradients end -# function MxVDerivative!(Phi_prime::Vector{Float64}, Phi::Vector{Float64}, -# D::Matrix{Float64}, N::Int) -# for i in 0:N -# t = 0.0 -# for j in 0:N -# t += D[i + 1, j + 1] * Phi[j + 1] # Adjust for 1-based indexing -# end -# Phi_prime[i + 1] = t -# end -# end - function compute_flux_array!(Flux, u, direction, equations) Np, nvars = size(Flux) for i in 1:Np From 44c3f5fb0baeb3e1b2d66e5264ef60c095a2571b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 14 Jul 2025 10:59:40 +0100 Subject: [PATCH 34/38] Use components of derivative vector for computing the vorticity. --- src/callbacks_step/analysis_dg2d.jl | 6 ++---- src/equations/compressible_euler_multicomponent_2d.jl | 4 ++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 202ab5c7ab8..48bcdfdae7e 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -260,15 +260,13 @@ function integrate(func::Func, u, if (m[1].nargs == 4) && (func != cons2cons) gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations) - gradients = Vector([gradients_x, gradients_y]) - @autoinfiltrate return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) + gradients_local = Vector([gradients_x[:, i, j, element], gradients_y[:, i, j, element]]) - func(u_local, gradients, equations) - return 0.0 + func(u_local, gradients_local, equations) end end end diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 844ea6e4e33..3167442a42c 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -878,8 +878,8 @@ end @inline function vorticity(u, gradients, equations::CompressibleEulerMulticomponentEquations2D) # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. - _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) - _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) + dv1dx, dv2dx = gradients[1][1:2] + dv1dy, dv2dy = gradients[2][1:2] return dv2dx - dv1dy end From d825c9313a3618f10259e3013bf3a98121b72286 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 14 Jul 2025 11:31:09 +0100 Subject: [PATCH 35/38] Removed factor rho from the enstrophy. --- src/equations/compressible_euler_multicomponent_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 3167442a42c..2c3cff9fe1b 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -872,7 +872,7 @@ end omega = vorticity(u, gradients, equations) - return 0.5f0 * u[1] * omega^2 + return 0.5f0 * omega^2 end @inline function vorticity(u, gradients, From 3cdbc821a8ae18c748feffcebb066d78b9977bc7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 24 Jul 2025 11:27:16 +0100 Subject: [PATCH 36/38] Replaced DG derivative with polynoial interpolation. --- src/callbacks_step/analysis_dg2d.jl | 94 +++++++++++++++++++++++++---- 1 file changed, 82 insertions(+), 12 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index dea699e7c66..af803fed9f2 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -258,8 +258,8 @@ function integrate(func::Func, u, end end if (m[1].nargs == 4) && (func != cons2cons) - gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations) - gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations) + gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations, mesh) + gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations, mesh) return integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg @@ -271,30 +271,100 @@ function integrate(func::Func, u, end end +# Return the derivatives of the interpolatied polynomial. +function lagrange_derivatives(x, y) + n = length(x) + dydx = zeros(n) + + for i in 1:n + for j in 1:n + if i != j + term = y[j] / (x[i] - x[j]) + for k in 1:n + if k != i && k != j + term *= (x[i] - x[k]) / (x[j] - x[k]) + end + end + dydx[i] += term + end + end + end + + return dydx +end + # Andrew's functions for computing the derivatives. function DGSpaceDerivative_WeakForm!(dg, cache, u, direction::Int, - equations) + equations, + mesh) # Get the required variables. - @unpack derivative_dhat = dg.basis + @unpack derivative_dhat, weights = dg.basis + + # Compute the surface flux terms. + surface_flux_values = similar(u) + Trixi.calc_surface_integral!(surface_flux_values, u, mesh, equations, dg.surface_integral, dg, cache) # Translations. D = derivative_dhat spA_Dhat = D - J = 1 ./ cache.elements.inverse_jacobian - weights = dg.basis.weights - l_minus = dg.basis.boundary_interpolation[:, 1] - l_plus = dg.basis.boundary_interpolation[:, 2] n_vars, Np, _, n_elements = size(u) gradients = similar(u) - for elem in 1:n_elements - for v in 1:n_vars - # Contract D in the ξ̂ (first spatial) direction - gradients[v, :, :, elem] .= D * u[v, :, :, elem] + u_local = zeros((n_vars, Np, Np)) + + if direction == 1 + for elem in 1:n_elements + # We work in primitive variables here. + for i in 1:Np + for j in 1:Np + u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations) + end + end + # Volume contributions (weak form) + for v in 1:n_vars + # Interpolate polynomial. + for j in 1:Np + x = cache.elements.node_coordinates[1, :, j, elem] + y = u_local[v, :, j] + gradients[v, :, j, elem] .= lagrange_derivatives(x, y) + end + +# # Contract D in the ξ̂ (first spatial) direction +# gradients[v, :, :, elem] .= D * u[v, :, :, elem] +# # Average over element. + + average = sum(gradients[v, :, :, elem]/Np^2) + gradients[v, :, :, elem] .= average + end + end + elseif direction == 2 + for elem in 1:n_elements + # We work in primitive variables here. + for i in 1:Np + for j in 1:Np + u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations) + end + end + # Volume contributions (weak form) + for v in 1:n_vars + # Interpolate polynomial. + for i in 1:Np + x = cache.elements.node_coordinates[2, i, :, elem] + y = u_local[v, i, :] + gradients[v, i, :, elem] .= lagrange_derivatives(x, y) + end + +# # Contract D in the ξ̂ (first spatial) direction +# gradients[v, :, :, elem] .= u[v, :, :, elem] * D' + + # Average over element. + average = sum(gradients[v, :, :, elem]/Np^2) + gradients[v, :, :, elem] .= average + end end end From 52b267b6453e6bd9c6d487ef03d334acf5768483 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 10 Oct 2025 15:11:26 +0100 Subject: [PATCH 37/38] Added max_abs_speed_naive for CompressibleEulerMulticomponentEquations2D and normal_direction::AbstractVector. --- .../compressible_euler_multicomponent_2d.jl | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index cf8299ef2dd..df088ce304a 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -681,6 +681,44 @@ end return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + + # Normalize the normal vector (if not guaranteed to be unit) + n1, n2 = normal_direction + norm_n = sqrt(n1^2 + n2^2) + n1 /= norm_n + n2 /= norm_n + + # Densities and gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Velocities + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + # Normal velocity components + v_n_ll = v1_ll * n1 + v2_ll * n2 + v_n_rr = v1_rr * n1 + v2_rr * n2 + + # Pressure and sound speeds + p_ll = (gamma_ll - 1f0) * (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) + p_rr = (gamma_rr - 1f0) * (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) + + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + return max(abs(v_n_ll), abs(v_n_rr)) + max(c_ll, c_rr) +end + # Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` @inline function max_abs_speed(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) From f891791a3e42008754f587b84de1f9a0fc6733fb Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 4 Nov 2025 16:08:39 +0000 Subject: [PATCH 38/38] Applied auto formatter. --- .../compressible_euler_multicomponent_2d.jl | 75 +++++++++++-------- 1 file changed, 45 insertions(+), 30 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index df088ce304a..156aeda2e48 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -61,12 +61,10 @@ struct CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real} < function CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}(gammas::SVector{NCOMP, RealT}, gas_constants::SVector{NCOMP, - RealT}) where { - NVARS, + RealT}) where {NVARS, NCOMP, RealT <: - Real - } + Real} NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `gas_constants` have to be filled with at least one value")) @@ -93,11 +91,9 @@ function CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants) end @inline function Base.real(::CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, - RealT}) where { - NVARS, + RealT}) where {NVARS, NCOMP, - RealT - } + RealT} RealT end @@ -143,7 +139,8 @@ function initial_condition_convergence_test(x, t, prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * rho / (1 - 2^ncomponents(equations)) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) prim1 = rho * v1 prim2 = rho * v2 @@ -186,7 +183,8 @@ Source terms used for convergence tests in combination with du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * tmp1 / (1 - 2^ncomponents(equations)) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) du1 = tmp5 du2 = tmp5 @@ -225,7 +223,8 @@ function initial_condition_weak_blast_wave(x, t, RealT(1.1691) / (1 - 2^ncomponents(equations)) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi @@ -418,10 +417,12 @@ Adaption of the entropy conserving two-point flux by rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) # Iterating over all partial densities rho_ll = density(u_ll, equations) @@ -459,7 +460,8 @@ Adaption of the entropy conserving two-point flux by help2 = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) for i in eachcomponent(equations) help1 += f_rho[i] * cv[i] help2 += f_rho[i] @@ -470,7 +472,8 @@ Adaption of the entropy conserving two-point flux by v2_avg * f2 else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) for i in eachcomponent(equations) help1 += f_rho[i] * cv[i] help2 += f_rho[i] @@ -508,10 +511,12 @@ See also rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) # Iterating over all partial densities rho_ll = density(u_ll, equations) @@ -554,7 +559,8 @@ See also f_rho_sum = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -564,7 +570,8 @@ See also 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -588,10 +595,12 @@ end rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) # Iterating over all partial densities rho_ll = density(u_ll, equations) @@ -636,7 +645,8 @@ end f_rho_sum = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -710,8 +720,8 @@ end v_n_rr = v1_rr * n1 + v2_rr * n2 # Pressure and sound speeds - p_ll = (gamma_ll - 1f0) * (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) - p_rr = (gamma_rr - 1f0) * (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) + p_ll = (gamma_ll - 1.0f0) * (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) + p_rr = (gamma_rr - 1.0f0) * (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) c_ll = sqrt(gamma_ll * p_ll / rho_ll) c_rr = sqrt(gamma_rr * p_rr / rho_rr) @@ -813,7 +823,8 @@ end rho_v1, rho_v2, rho_e = u prim_rho = SVector{ncomponents(equations), real(equations)}(u[i + 3] - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) rho = density(u, equations) v1 = rho_v1 / rho @@ -856,7 +867,8 @@ end gas_constants[i] * (1 + log(u[i + 3])) - v_square / (2 * T)) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) w1 = gas_constant * v1 * rho_p w2 = gas_constant * v2 * rho_p @@ -881,7 +893,8 @@ end (2 * T)) / gas_constants[i] - 1) - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) RealT = eltype(w) rho = zero(RealT) @@ -908,7 +921,8 @@ end v1, v2, p = prim cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i + 3] - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) rho = density(prim, equations) gamma = totalgamma(prim, equations) @@ -1001,7 +1015,8 @@ end @inline function densities(u, v, equations::CompressibleEulerMulticomponentEquations2D) return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v - for i in eachcomponent(equations)) + for i in + eachcomponent(equations)) end @inline function velocity(u, equations::CompressibleEulerMulticomponentEquations2D) @@ -1019,7 +1034,7 @@ end end @inline function enstrophy_multi_euler(u, gradients, - equations::CompressibleEulerMulticomponentEquations2D) + equations::CompressibleEulerMulticomponentEquations2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations)