added slurm job id to file output, fixed y_save and x_save scope issue.

plumeSSME.jl

@@ -41,6 +41,8 @@ s = parse(Int, ARGS[2]) #s steps in y direction
 ψ_mult= parse(Float64, ARGS[5])
 ϕ_mult= parse(Float64, ARGS[6])
 job_id = ARGS[7]
+set_T = ARGS[8] #imported or ambient
+set_u = ARGS[9] #imported or ambient
 
 println("Reading arguments from bash: s = ", s, " n = ", n)
 println("Reading arguments from bash: ψ_init = ", ψ_init, " ϕ_init = ", ϕ_init)
@@ -49,20 +51,21 @@ println("Reading arguments from bash: ψ_mult = ", ψ_mult, " ϕ_mult = ", ϕ_mu
 T_save = zeros(length(h), s, n)
 u_save = zeros(length(h), s, n)
 χ_save = zeros(length(h), s, n, n_species)
-x_save = zeros(n)
-y_save = zeros(s)
+y = zeros(s,n)
+x = zeros(n)
+
 gas_g = StructArray{gas_type}(undef,s,n,length(h))
 gas_g .= [gas_type(0)]
 
 println("started populating phi and psi")
 Δϕ = ϕ_init * ones(n) #step size in phi
 for i = 2:n
-    Δϕ[i] = ϕ_mult * Δϕ[i-1] #enlarge with each step by 1.1
+    Δϕ[i] = ϕ_mult * Δϕ[i-1] #enlarge with each step by mult
 end
 
 Δψ = ψ_init * ones(s) #s vertical grid points
 for i = 2:s
-    Δψ[i] = ψ_mult * Δψ[i-1] #enlarge with each step by 1.1
+    Δψ[i] = ψ_mult * Δψ[i-1] #enlarge with each step by mult
 end
 
 println("started altitude for loop")
@@ -71,21 +74,27 @@ for m = 1:lastindex(h)
 
     ### IMPORT SHOCK EXIT CONDITIONS ###
 
-    u0 = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/u")#initial plume velocity
-    u0 = convert(AbstractFloat, u0[2])
-
-    T0 = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/T") #initial plume temperature
-    T0 = convert(AbstractFloat, T0[2])
-
     p_all = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/P")
     p = convert(AbstractFloat, p_all[2])
 
     u_a = 1.11849E-19 * big(h[m])^5 - 1.14814E-14 * big(h[m])^4 + 4.22542E-10 * big(h[m])^3 - 6.92322E-06 * big(h[m])^2 + 6.58761E-02 * big(h[m]) + 5.37920E+01
 
     #curve fit #a = ambient vel [m/s] (speed of rocket) 
     T_a = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/T_a")
-    #T0 = T_a #for testing see vincent messages
-    #u0 = u_a #for testing see vincent messages
+
+    if set_T == "ambient"
+        T0 = T_a #for testing
+    else
+        T0 = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/T") #initial plume temperature
+        T0 = convert(AbstractFloat, T0[2])
+    end
+
+    if set_u == "ambient"
+        u0 = u_a #for testing
+    else
+        u0 = HDF5.h5read("/home/chinahg/GCresearch/rocketemissions/plot_data.h5", h_string * "m/u")#initial plume velocity
+        u0 = convert(AbstractFloat, u0[2])
+    end
 
     ### CALCULATE VELOCITY AND TEMPERATURE FIELDS (NO CHEMISTRY) ###
 
@@ -98,16 +107,16 @@ for m = 1:lastindex(h)
     y_init = compute_y(u_init, T_init, Δψ, R, p)
 
     #Geometry of plume
-    radius = 1.147
+    radius = 0.1#1.147
     u_init[y_init.>radius] .= convert(AbstractFloat, u_a)
     T_init[y_init.>radius] .= convert(AbstractFloat, T_a)
 
     ambient = AmbientConditions(u_a, T_a, 1.0, 1.0, R, p, T0, u0)
 
     #Solve for T and u at all steps: NO CHEMISTRY
     x, y, u, T, ϵ = solve_exhaust_flow(u_init, T_init, ambient, n, Δϕ, Δψ)
-    x_save = x
-    y_save = y
+    println("solve exhaust: ", size(y))
+
     T_save[m, :, :] = T
     u_save[m, :, :] = u
 
@@ -177,7 +186,7 @@ for m = 1:lastindex(h)
     #dummy_reactor = ct.IdealGasReactor(gas)
 
     println("starting splitting")
-    Threads.@threads for i = 1:n-1 #x
+    for i = 1:n-1 #x
 
         for j = 1:n_species #species
             #calculate f0 at half step 0.5*Δϕ (x)
@@ -216,19 +225,18 @@ if plot_flag
         T = T_save[m, :, :] #temperature
         u = u_save[m, :, :] #velocity
         χ = χ_save[m, :, :, :] #concentrations
-        x = x_save
-        y = y_save
+        println("before regrid: ", size(y))
         h_string = string(h[m])
         P_atm = 101325 #placeholder, need altitude dependent
         R = 287 #[J/kgK] placeholder
 
         #REGRID SOLUTION
-        #xx, yy, u_g, T_g, χ_gO2 =  regrid_solution(x, y, u, T, χ[:,:,4], 0.01)
-        #xx, yy, u_g, T_g, χ_gN2 =  regrid_solution(x, y, u, T, χ[:,:,48], 0.01)
+        xx, yy, u_g, T_g, χ_gO2 =  regrid_solution(x, y, u, T, χ[:,:,4], 0.01)
+        xx, yy, u_g, T_g, χ_gN2 =  regrid_solution(x, y, u, T, χ[:,:,48], 0.01)
         #xx, yy, u_g, T_g, χ_gNO2 = regrid_solution(x, y, u, T, χ[:, :, 37], 0.01)
         #xx, yy, u_g, T_g, χ_gN2O = regrid_solution(x, y, u, T, χ[:, :, 38], 0.01)
         #xx, yy, u_g, T_g, χ_gNO = regrid_solution(x, y, u, T, χ[:, :, 36], 0.01)
-        #xx, yy, u_g, T_g, χ_gAr = regrid_solution(x, y, u, T, χ[:, :, 49], 0.01)
+        xx, yy, u_g, T_g, χ_gAr = regrid_solution(x, y, u, T, χ[:, :, 49], 0.01)
 
         #PLOT 2D MAPS AND SAVE
         x_max = 20 #[m]
@@ -304,22 +312,22 @@ if plot_flag
         mdot_N2[s,:] = mdot_N2[s-1,:]
         mdot_O2[s,:] = mdot_O2[s-1,:]
 
-        #fig,axc = plt.subplots()
-        #axc.plot(ϕ, u[s,:], label="5000/5000")
-        #axc.plot(ϕ, u[s-1,:], label="4999/5000")
-        #axc.plot(ϕ, u[s-50,:], label="4950/5000")
-        #axc.plot(ϕ, u[s-100,:], label="4900/5000")
-        #axc.plot(ϕ, u[s-150,:], label="4850/5000")
-        #axc.plot(ϕ, u[s-200,:], label="4800/5000")
-        #axc.plot(ϕ, u[s-250,:], label="4750/5000")
-        #axc.plot(ϕ, u[s-1000,:], label="4000/5000")
+        # fig,axc = plt.subplots()
+        # axc.plot(ϕ, u[s,:], label="5000/5000")
+        # axc.plot(ϕ, u[s-1,:], label="4999/5000")
+        # axc.plot(ϕ, u[s-50,:], label="4950/5000")
+        # axc.plot(ϕ, u[s-100,:], label="4900/5000")
+        # axc.plot(ϕ, u[s-150,:], label="4850/5000")
+        # axc.plot(ϕ, u[s-200,:], label="4800/5000")
+        # axc.plot(ϕ, u[s-250,:], label="4750/5000")
+        # axc.plot(ϕ, u[s-1000,:], label="4000/5000")
 
-        #axc.set_xlabel("ϕ")
-        #axc.set_ylabel("m/s at ψ = 0.1")
-        #axc.set_xscale("log")
-        #legend = axc.legend(loc = "upper right")
-        #fig.suptitle("Centerline Velocity")
-        #savefig("/home/chinahg/GCresearch/rocketemissions/rockettests/" * h_string * "m/" * job_id * "_u_cent.png")
+        # axc.set_xlabel("ϕ")
+        # axc.set_ylabel("m/s at ψ = 0.1")
+        # axc.set_xscale("log")
+        # legend = axc.legend(loc = "upper right")
+        # fig.suptitle("Centerline Velocity")
+        # savefig("/home/chinahg/GCresearch/rocketemissions/rockettests/" * h_string * "m/" * job_id * "_u_cent.png")
 
         fig,axR = plt.subplots(2,2)
         axR[1,1].plot(ϕ, mdot_Ar[1,:])
@@ -357,33 +365,33 @@ if plot_flag
         savefig("/home/chinahg/GCresearch/rocketemissions/rockettests/" * h_string * "m/" * job_id * "_u.png")
 
 
-        #fig,axs = plt.subplots(2,2)
-        #axs[1,1].plot(ϕ, mdot_tot)
-        #axs[1,1].set_xlabel("ϕ")
-        #axs[1,1].set_ylabel("Total mass flow rate [kg/s]")
-        #axs[1,1].set_title("Total mdot")
-        #axs[1,1].set_xscale("log")
-    #
-        #axs[1,2].plot(ϕ, mdot_N2_sumTrunc)
-        #axs[1,2].set_xlabel("ϕ")
-        #axs[1,2].set_ylabel("N2 mass flow rate [kg/s]")
-        #axs[1,2].set_title("N2")
-        #axs[1,2].set_xscale("log")
-    #
-        #axs[2,1].plot(ϕ, mdot_O2_sumTrunc)
-        #axs[2,1].set_xlabel("ϕ")
-        #axs[2,1].set_ylabel("O2 mass flow rate [kg/s]")
-        #axs[2,1].set_xscale("log")
-    #
-        #axs[2,2].plot(ϕ, mdot_Ar_sumTrunc)
-        #axs[2,2].set_xlabel("ϕ")
-        #axs[2,2].set_ylabel("Ar mass flow rate [kg/s]")
-        #axs[2,2].set_title("Ar")
-        #axs[2,2].set_xscale("log")
-        #
-        #fig.suptitle("Truncated Mass Flow Rates, s = 250/300")
-        #fig.tight_layout()
-        #savefig("/home/chinahg/GCresearch/rocketemissions/rockettests/" * h_string * "m/" * job_id * "_mdot_trunc.png")
+        # fig,axs = plt.subplots(2,2)
+        # axs[1,1].plot(ϕ, mdot_tot)
+        # axs[1,1].set_xlabel("ϕ")
+        # axs[1,1].set_ylabel("Total mass flow rate [kg/s]")
+        # axs[1,1].set_title("Total mdot")
+        # axs[1,1].set_xscale("log")
+
+        # axs[1,2].plot(ϕ, mdot_N2_sumTrunc)
+        # axs[1,2].set_xlabel("ϕ")
+        # axs[1,2].set_ylabel("N2 mass flow rate [kg/s]")
+        # axs[1,2].set_title("N2")
+        # axs[1,2].set_xscale("log")
+
+        # axs[2,1].plot(ϕ, mdot_O2_sumTrunc)
+        # axs[2,1].set_xlabel("ϕ")
+        # axs[2,1].set_ylabel("O2 mass flow rate [kg/s]")
+        # axs[2,1].set_xscale("log")
+
+        # axs[2,2].plot(ϕ, mdot_Ar_sumTrunc)
+        # axs[2,2].set_xlabel("ϕ")
+        # axs[2,2].set_ylabel("Ar mass flow rate [kg/s]")
+        # axs[2,2].set_title("Ar")
+        # axs[2,2].set_xscale("log")
+
+        # fig.suptitle("Truncated Mass Flow Rates, s = 250/300")
+        # fig.tight_layout()
+        # savefig("/home/chinahg/GCresearch/rocketemissions/rockettests/" * h_string * "m/" * job_id * "_mdot_trunc.png")
 
         fig,axX = plt.subplots(2,2)
         axX[1,1].plot(ϕ, mdot_tot)

plumefunctions.jl

@@ -138,9 +138,12 @@ grid spacing in ϕ-ψ space.
 """
 function compute_y(u::AbstractVector, T::AbstractVector, Δψ::AbstractVector,
     R::AbstractFloat, p::AbstractFloat)
-
+    try 
+        (2 * cumsum((R / p) * T ./ u) .* Δψ) .^ 0.5
+    catch
+        println(T,u)
+    end
     return (2 * cumsum((R / p) * T ./ u) .* Δψ) .^ 0.5
-
 end
 
 """
@@ -243,9 +246,10 @@ Maps back solution from a grid in ϕ-ψ space to a grid in x-y space.
 - `T` Temperature values of solution
 - `y_spacing` Desired spacing in `y` direction for output grid
 """
-function regrid_solution(x::AbstractVector, y::AbstractMatrix, u::AbstractMatrix,
-    T::AbstractMatrix, χ::AbstractMatrix, y_spacing::AbstractFloat)
+function regrid_solution(x::Array, y::Array, u::Array,
+    T::Array, χ::Array, y_spacing::Float64)
 
+    println("regrid: ",size(y))
     yy = 0:y_spacing:maximum(y)
 
     u_gridded = zeros((size(yy)[1], size(x)[1]))

slurm/plume.sh

@@ -2,36 +2,41 @@
 
 
 #SBATCH --time=48:00:00
-#SBATCH --job-name="2000x5000"
+#SBATCH --job-name="1000x2000"
 #SBATCH [email protected]
 #SBATCH --mail-type=BEGIN,END
 #SBATCH -o slurm-%j.out
 #xSBATCH -e slurm-%j.err
 #SBATCH --nodes=1
 #SBATCH --ntasks-per-node=1
-#SBATCH --cpus-per-task=8
+#SBATCH --cpus-per-task=1
 #SBATCH --partition=normal
-#SBATCH --mem=56000MB
+#SBATCH --mem=50000MB
 #####################################
 
-n=$11000 #2000
-s=$25000 #5000
-psi_init=$30.001 #0.001
-phi_init=$40.001 #0.01
-psi_mult=$51.01 #1.01
-phi_mult=$61.01 #1.01
-job_id=$7#SLURM_JOB_ID
+n=$15000
+s=$25000
+psi_init=$30.001
+phi_init=$40.001
+psi_mult=$51.01
+phi_mult=$61.01
+job_id=$SLURM_JOBID
+set_T=$8imported #imported or ambient
+set_u=$9imported #imported or ambient
+
+echo "Job ID: $job_id"
 echo "Setting s = $s and n = $n"
 echo "psi_init = $psi_init, phi_init = $phi_init"
 echo "psi_mult = $psi_mult, phi_init = $phi_mult"
+echo "Temperature = $set_T, Velocity = $set_u"
 
 now=$(date +"%T")
 echo "Start time : $now"
 
 echo "Number of CPUs per task: $SLURM_CPUS_PER_TASK"
-julia --project=/home/chinahg/GCresearch/rocketemissions/ --threads $SLURM_CPUS_PER_TASK /home/chinahg/GCresearch/rocketemissions/plumeSSME.jl $n $s $psi_init $phi_init $psi_mult $phi_mult $job_id
+julia --project=/home/chinahg/GCresearch/rocketemissions/  /home/chinahg/GCresearch/rocketemissions/plumeSSME.jl $n $s $psi_init $phi_init $psi_mult $phi_mult $job_id $set_T $set_u
 
 now=$(date +"%T")
 echo "End time : $now"
 
-#
+#--threads $SLURM_CPUS_PER_TASK