testing the dependency on dt, drawing frequency

docs/src/API.md

@@ -93,6 +93,7 @@ HetIceNucleation.ABIFM_J
 HetIceNucleation.P3_deposition_N_i
 HetIceNucleation.P3_het_N_i
 HetIceNucleation.INP_concentration_frequency
+HetIceNucleation.INP_concentration_mean
 ```
 
 # Homogeneous ice nucleation

parcel/Frostenberg_different_dt.jl

@@ -0,0 +1,112 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+FT = Float32
+# get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+
+# Initial conditions
+ρₗ = wps.ρw
+Nₐ = FT(0)
+Nₗ = FT(500 * 1e3)
+Nᵢ = FT(0)
+r₀ = FT(1e-6)
+p₀ = FT(800 * 1e2)
+T₀ = FT(251)
+qᵥ = FT(8.1e-4)
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
+qᵢ = FT(0)
+x_sulph = FT(0.01)
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_v = TD.Parameters.R_v(tps)
+Rₐ = TD.gas_constant_air(tps, q)
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+e = eᵥ(qᵥ, p₀, Rₐ, R_v)
+Sₗ = FT(e / eₛ)
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+# Simulation parameters passed into ODE solver
+w = FT(0.7)                                # updraft speed
+t_max = FT(1200)
+aerosol = CMP.Illite(FT)
+heterogeneous = "Frostenberg_random"
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+DSD = "Monodisperse"
+const_dt_range = range(FT(0.1), stop = FT(1), length = 4) #changing the time step and drawing in every time step
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]")
+ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice")
+
+for const_dt in const_dt_range
+
+    params = parcel_params{FT}(
+        const_dt = const_dt,
+        w = w,
+        aerosol = aerosol,
+        heterogeneous = heterogeneous,
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        size_distribution = DSD,
+        drawing_interval = const_dt,    #drawing in every time step
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, params)
+
+    # Plot results
+    MK.lines!(
+        ax1,
+        sol.t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        label = "dt=" * string(const_dt) * " s",
+    )
+    MK.lines!(
+        ax2,
+        sol.t / 60,
+        sol[3, :],
+        label = "dt=" * string(const_dt) * " s",
+    )
+    MK.lines!(
+        ax3,
+        sol.t / 60,
+        sol[6, :] * 1e3,
+        label = "dt=" * string(const_dt) * " s",
+    )
+    MK.lines!(
+        ax4,
+        sol.t / 60,
+        sol[5, :] * 1e3,
+        label = "dt=" * string(const_dt) * " s",
+    )
+    sol_Nₗ = sol[8, :]
+    sol_Nᵢ = sol[9, :]
+    MK.lines!(ax5, sol.t / 60, sol_Nₗ, label = "dt=" * string(const_dt) * " s")
+    MK.lines!(ax6, sol.t / 60, sol_Nᵢ, label = "dt=" * string(const_dt) * " s")
+end
+
+for ax in [ax1, ax2, ax3, ax4, ax5, ax6]
+    MK.axislegend(
+        ax,
+        framevisible = false,
+        labelsize = 12,
+        orientation = :horizontal,
+        nbanks = 2,
+        position = :lc,
+    )
+end
+
+MK.save("dt.svg", fig)

parcel/Frostenberg_different_frequency.jl

@@ -0,0 +1,135 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+FT = Float32
+# get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+
+# Initial conditions
+ρₗ = wps.ρw
+Nₐ = FT(0)
+Nₗ = FT(500 * 1e3)
+Nᵢ = FT(0)
+r₀ = FT(1e-6)
+p₀ = FT(800 * 1e2)
+T₀ = FT(251)
+qᵥ = FT(8.1e-4)
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
+qᵢ = FT(0)
+x_sulph = FT(0.01)
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_v = TD.Parameters.R_v(tps)
+Rₐ = TD.gas_constant_air(tps, q)
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+e = eᵥ(qᵥ, p₀, Rₐ, R_v)
+Sₗ = FT(e / eₛ)
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+# Simulation parameters passed into ODE solver
+w = FT(0.7)                                # updraft speed
+t_max = FT(1200)
+aerosol = CMP.Illite(FT)
+heterogeneous = "Frostenberg_random"
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+DSD = "Monodisperse"
+const_dt = FT(1)                                                 #constent time step
+drawing_interval_range = range(FT(1), stop = FT(5), length = 4)  # changing the time between drawing
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]")
+ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice")
+
+for drawing_interval in drawing_interval_range
+
+    params = parcel_params{FT}(
+        const_dt = const_dt,
+        w = w,
+        aerosol = aerosol,
+        heterogeneous = heterogeneous,
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        size_distribution = DSD,
+        drawing_interval = drawing_interval,
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, params)
+
+    # Plot results
+    MK.lines!(
+        ax1,
+        sol.t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+    MK.lines!(
+        ax2,
+        sol.t / 60,
+        sol[3, :],
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+    MK.lines!(
+        ax3,
+        sol.t / 60,
+        sol[6, :] * 1e3,
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+    MK.lines!(
+        ax4,
+        sol.t / 60,
+        sol[5, :] * 1e3,
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+    sol_Nₗ = sol[8, :]
+    sol_Nᵢ = sol[9, :]
+    MK.lines!(
+        ax5,
+        sol.t / 60,
+        sol_Nₗ,
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+    MK.lines!(
+        ax6,
+        sol.t / 60,
+        sol_Nᵢ,
+        label = "f =" *
+                string(round(const_dt / drawing_interval, digits = 1)) *
+                "/s",
+    )
+end
+
+for ax in [ax1, ax2, ax3, ax4, ax5, ax6]
+    MK.axislegend(
+        ax,
+        framevisible = false,
+        labelsize = 12,
+        orientation = :horizontal,
+        nbanks = 2,
+        position = :lc,
+    )
+end
+
+MK.save("frequency.svg", fig)