From 41c03059b733a9aeb70352bf0428481467ba8bf3 Mon Sep 17 00:00:00 2001 From: smarras79 Date: Thu, 18 Jul 2024 09:59:56 +0200 Subject: [PATCH] merged with Hang fixes to soundSpeed and cfl calculator for GPU --- .../equations/CompEuler/theta/user_inputs.jl | 1 - src/kernel/physics/constitutiveLaw.jl | 5 ++- src/kernel/physics/soundSpeed.jl | 43 +++++++++++-------- src/kernel/solvers/TimeIntegrators.jl | 12 +++--- 4 files changed, 36 insertions(+), 25 deletions(-) diff --git a/problems/equations/CompEuler/theta/user_inputs.jl b/problems/equations/CompEuler/theta/user_inputs.jl index 168a2bf1..badb7c6a 100644 --- a/problems/equations/CompEuler/theta/user_inputs.jl +++ b/problems/equations/CompEuler/theta/user_inputs.jl @@ -13,7 +13,6 @@ function user_inputs() #:tend => 1000.0, #:lrestart => true, :restart_input_file_path => "./output/CompEuler/theta/output-19Nov2023-115126", - #:ndiagnostics_outputs => 2, :diagnostics_at_times => (100, 200, 300, 400, 500, 600, 700, 800, 900, 1000), :case => "rtb", :lsource => true, diff --git a/src/kernel/physics/constitutiveLaw.jl b/src/kernel/physics/constitutiveLaw.jl index ad354781..3c6dbe3b 100644 --- a/src/kernel/physics/constitutiveLaw.jl +++ b/src/kernel/physics/constitutiveLaw.jl @@ -17,8 +17,11 @@ end function perfectGasLaw_ρθtoP(PhysConst::PhysicalConst; ρ=1.25, θ=300.0) return PhysConst.C0*(ρ*θ)^PhysConst.γ #Press - #return PhysConst.pref*(ρ*θ*PhysConst.Rair/PhysConst.pref)^PhysConst.cpovercv #Press + +end +function perfectGasLaw_ρθtoP(PhysConst::PhysicalConst, ρ::AbstractArray, θ::AbstractArray) + return PhysConst.C0 .* (ρ .* θ) .^ PhysConst.γ end function perfectGasLaw_ρθtoP!(Press::Array{Float64}, PhysConst::PhysicalConst; ρ=1.25, θ=300.0) diff --git a/src/kernel/physics/soundSpeed.jl b/src/kernel/physics/soundSpeed.jl index 0ecd330c..abfb744f 100644 --- a/src/kernel/physics/soundSpeed.jl +++ b/src/kernel/physics/soundSpeed.jl @@ -1,20 +1,27 @@ -function soundSpeed(npoin, integrator) - - #speed of sound - PhysConst = PhysicalConst{TFloat}() - ctmp = Float32(0.0) - c = Float32(0.0) - for i=1:npoin - ρ = integrator.u[i] - θ = integrator.u[3*npoin+i] - p = perfectGasLaw_ρθtoP(PhysConst; ρ=ρ, θ=θ) - c = sqrt(PhysConst.γ*p/ρ) - c = max(c, ctmp) - ctmp = c - end +function soundSpeed(npoin, integrator, SD) - return c + # Physical constants + PhysConst = PhysicalConst{Float32}() + pos::TInt = 2 + if (SD == NSD_2D()) + pos = 3 + elseif (SD == NSD_3D()) + pos = 4 + end + # Initialize arrays + ρ = integrator.u[1:npoin] + θ = integrator.u[pos*npoin+1:(pos+1)*npoin] + # Compute pressure using vectorized operation + p = perfectGasLaw_ρθtoP(PhysConst, ρ, θ) + + # Compute speed of sound using vectorized operation + c = sqrt.(PhysConst.γ .* p ./ ρ) + + # Find the maximum speed of sound + max_c = maximum(c) + + return max_c end function computeCFL(npoin, dt, Δs, integrator, SD::NSD_2D) @@ -32,7 +39,7 @@ function computeCFL(npoin, dt, Δs, integrator, SD::NSD_2D) velomax = max(umax, vmax) #speed of sound - c = soundSpeed(npoin, integrator) + c = soundSpeed(npoin, integrator, SD) cfl_u = velomax*dt/Δs #Advective CFL cfl_c = c*dt/Δs #Acoustic CFL @@ -55,13 +62,13 @@ function computeCFL(npoin, dt, Δs, integrator, SD::NSD_3D) #w ieq = 4 idx = (ieq-1)*npoin - vmax = maximum(integrator.u[idx+1:ieq*npoin]) + wmax = maximum(integrator.u[idx+1:ieq*npoin]) #velomax velomax = max(umax, vmax, wmax) #speed of sound - c = soundSpeed(npoin, integrator) + c = soundSpeed(npoin, integrator, SD) cfl_u = velomax*dt/Δs #Advective CFL cfl_c = c*dt/Δs #Acoustic CFL diff --git a/src/kernel/solvers/TimeIntegrators.jl b/src/kernel/solvers/TimeIntegrators.jl index 9077a17d..3ddec26a 100644 --- a/src/kernel/solvers/TimeIntegrators.jl +++ b/src/kernel/solvers/TimeIntegrators.jl @@ -10,7 +10,7 @@ function time_loop!(inputs, params, u) params); #------------------------------------------------------------------------ - # Callback to plot on the run + # Runtime callbacks #------------------------------------------------------------------------ dosetimes = inputs[:diagnostics_at_times] idx_ref = Ref{Int}(0) @@ -29,8 +29,11 @@ function time_loop!(inputs, params, u) idx = idx_ref[] println(" # t=", integrator.t) + + #CFL computeCFL(params.mesh.npoin, inputs[:Δt], params.mesh.Δeffective_s, integrator, params.SD) - + + #Write results to file write_output(params.SD, integrator.u, integrator.t, idx, params.mesh, inputs[:output_dir], inputs, @@ -40,10 +43,9 @@ function time_loop!(inputs, params, u) end - cb = DiscreteCallback(condition, affect!) - + cb = DiscreteCallback(condition, affect!) #------------------------------------------------------------------------ - # END Callback to plot on the run + # END runtime callbacks #------------------------------------------------------------------------ @time solution = solve(prob,