merged with Hang fixes to soundSpeed and cfl calculator for GPU

smarras79 · Jul 18, 2024 · 41c0305 · 41c0305
1 parent 61ff1c0
commit 41c0305
Show file tree

Hide file tree

Showing 4 changed files with 36 additions and 25 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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,