docs suggestions and cleanup in stochastic parcel

.gitignore

@@ -21,6 +21,11 @@
 !docs/src/assets/*.svg
 docs/build/
 docs/site/
+# ignore output from Literate.jl
+docs/src/guides/literated/
 
 # Deps
 Manifest.toml
+
+# Settings
+.vscode/

docs/Project.toml

@@ -12,6 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 

docs/src/IceNucleationParcel0D.md

@@ -353,46 +353,56 @@ Three different implementations of this parametrization are used in the parcel m
   with the mean and standard deviation defined in [Frostenberg2023](@cite).
   The inverse timescale of the process is set by ``\gamma``.
 
-The stochastic implementation is based on the equation for a generic
-  stochastic process ``x`` with Gaussian random noise (i.e. Wiener process):
+The stochastic implementation is based on the [Ornstein-Uhlenbeck process](https://en.wikipedia.org/wiki/Ornstein–Uhlenbeck_process), 
+in which the variable ``x`` is a mean-reverting process perturbed by Gaussian random noise (i.e. increments of the Wiener process ``W``):
 ```math
 \begin{equation}
-  dx = - \gamma \; x \; dt + g dW; \;\;\;\;\;\;\;   dW = {\bf{N}}(0, dt^2),
+  dx = - \gamma(x - \mu)dt + \sqrt{2\gamma} \sigma dW; \quad\quad   dW \sim N(0, dt),
 \end{equation}
 ```
-where ``\bf{N}`` is the normal distribution.
-For constant ``\gamma`` and ``g``, ``x`` has solution
+where ``N`` is a zero-mean normal distribution with variance ``dt``.
+For constant ``\gamma`` and ``\sigma``, and given some initial condition ``x(0)=x_0``, ``x`` has the analytical solution:
 ```math
 \begin{equation}
-  x(t) = x_0 \; e^{-\gamma t} + g \; \int_0^t e^{\gamma(s-t)} dW
+  x(t) = 
+      x_0 e^{-\gamma t} + \mu (1 - e^{-\gamma t})
+    + \sqrt{2\gamma} \sigma \int_0^t e^{-\gamma(t-s)} dW,
 \end{equation}
 ```
-where ``1/\gamma`` is the assumed timescale of the process.
-We can calculate the variance ``\bf{V}`` as,
+where ``\tau \equiv 1 / \gamma`` is the assumed timescale of the process. The process mean is ``x_0 e^{-\gamma t} + \mu (1 - e^{-\gamma t})``. We can calculate the variance ``\mathbb{V}(t)`` as,
 ```math
 \begin{equation}
-  {\bf{V}}(t) = g^2 \int_0^t e^{2\gamma(s-t)} ds = \frac{g^2}{2\gamma}(1 - e^{-2\gamma t}).
+  \mathbb{V}(t) 
+  = 2\gamma \sigma^2 \int_0^t e^{-2\gamma(t-s)} ds 
+  = \frac{g^2}{2\gamma} \left( 1 - e^{-2\gamma t} \right).
 \end{equation}
 ```
 
-We use this process to model ``x=\log(\text{INPC})``,
-  which tends toward the mean ``\mu(T)``.
-If we denote ``\tau = 1 / \gamma`` as the process timescale and
-  ``\sigma = g / \sqrt{2 / \gamma}`` as the process uncertainty,
-  then
+We use this process to model ``x=\log(\text{INPC})``, which tends toward a temperature-dependent mean value ``\mu(T)``.
+The equation for ``\log(\text{INPC})`` is then,
 ```math
 \begin{equation}
-  \frac{d\log(\text{INPC})}{dt} = - \frac{\log(\text{INPC}) - μ}{\tau} + \sigma \sqrt{\frac{2}{\tau \; dt}} \; {\bf{N}}(0, 1))
+  d\log(\text{INPC}) = 
+    - \frac{\log(\text{INPC}) - μ}{\tau} dt
+    + \sigma \sqrt{\frac{2}{\tau}} dW
 \end{equation}
 ```
 
-The following plot shows resuls of the parcel model with the `mean` (black line),
-  `random` (dotted lines) and `stochastic` (solid lines) parameterization options.
-We show results for two sampling intervals ``\Delta t`` (random),
-  two process time scales ``\tau`` (stochastic),
-  and two model time steps `dt`.
+This equation is currently implemented with the simple [Euler-Maruyama method](https://en.wikipedia.org/wiki/Euler–Maruyama_method), which is the stochastic analog of the forward Euler method for (deterministic) ordinary differential equations, so that
+```math
+\begin{equation}
+  \log(\text{INPC})_{t+dt} = \log(\text{INPC})_{t}
+    - \gamma\left(\log(\text{INPC})_t - μ(T_t)\right) dt
+    + \sigma \sqrt{2\gamma dt} z_t
+\end{equation}
+```
+where ``z_t \sim N(0,1)`` is a standard normal random variable.
+
+The following plot shows resuls of the parcel model with the `mean` (black line), `random` (dotted lines) and `stochastic` (solid lines) parameterization options.
+We show results for two sampling intervals ``\Delta t`` (random), two process time scales ``\tau`` (stochastic), and two model time steps `dt`.
 
 ```@example
-include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl")
+using Suppressor: @suppress # hide
+@suppress include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl") # hide
 ```
 ![](frostenberg_immersion_freezing.svg)

parcel/Example_Frostenberg_Immersion_Freezing.jl

@@ -137,18 +137,26 @@ end
 
 # Plotting
 fig = MK.Figure(size = (1200, 900))
-ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
-ax2 = MK.Axis(fig[1, 2], ylabel = "T [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_ice [1/cm3]")
-ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_liq [1/cm3]")
+plot_theme = MK.Theme(Axis = (; xgridvisible = false, ygridvisible = false))
+MK.set_theme!(plot_theme)
+ax1 = MK.Axis(fig[1, 1]; ylabel = "Ice Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2]; ylabel = "T [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_ice [1/cm3]")
+ax6 = MK.Axis(fig[3, 2]; xlabel = "Time [min]", ylabel = "N_liq [1/cm3]")
+MK.xlims!.(fig.content, 0, t_max / 60)
 ax7 = MK.Axis(
-    fig[4, 1:2],
+    fig[4, 1:2];
     xlabel = "T [K]",
     ylabel = "INPC [1/m3]",
     yscale = log10,
 )
+# top axis with time
+ax7_top = MK.Axis(fig[4, 1:2]; xaxisposition = :top, xlabel = "Time [min]")
+MK.hidespines!(ax7_top)
+MK.hideydecorations!(ax7_top)
+MK.xlims!(ax7_top, 0, t_max / 60)
 
 colors = [:red, :green, :orange, :limegreen]
 colors_mean = [:black, :gray]
@@ -194,5 +202,7 @@ for (sol, t, ll, cl) in zip(results_mean, time_mean, labels_mean, colors_mean)
         color = cl,
     )
 end
+MK.xlims!(ax7, extrema(results_mean[2][3, :]) |> reverse)
 fig[1:2, 3] = MK.Legend(fig, ax6, framevisible = false)
 MK.save("frostenberg_immersion_freezing.svg", fig)
+nothing

parcel/ParcelModel.jl

@@ -192,8 +192,7 @@ function run_parcel(IC, t_0, t_end, pp)
 
     FT = eltype(IC)
 
-    println(" ")
-    println("Size distribution: ", pp.size_distribution)
+    info = "\nSize distribution: $(pp.size_distribution)\n"
     if pp.size_distribution == "Monodisperse"
         distr = Monodisperse{FT}()
     elseif pp.size_distribution == "Gamma"
@@ -202,9 +201,9 @@ function run_parcel(IC, t_0, t_end, pp)
         throw("Unrecognized size distribution")
     end
 
-    println("Aerosol :", chop(string(typeof(pp.aerosol)), head = 29, tail = 9))
+    info *= "Aerosol: $(chop( string(typeof(pp.aerosol)), head = 29, tail = 9))\n"
 
-    println("Deposition: ", pp.deposition)
+    info *= "Deposition: $(pp.deposition)\n"
     if pp.deposition == "None"
         dep_params = Empty{FT}()
     elseif pp.deposition == "MohlerAF"
@@ -219,7 +218,7 @@ function run_parcel(IC, t_0, t_end, pp)
         throw("Unrecognized deposition mode")
     end
 
-    println("Heterogeneous: ", pp.heterogeneous)
+    info *= "Heterogeneous: $(pp.heterogeneous)\n"
     if pp.heterogeneous == "None"
         imm_params = Empty{FT}()
     elseif pp.heterogeneous == "ABIFM"
@@ -237,7 +236,7 @@ function run_parcel(IC, t_0, t_end, pp)
         throw("Unrecognized heterogeneous mode")
     end
 
-    println("Homogeneous: ", pp.homogeneous)
+    info *= "Homogeneous: $(pp.homogeneous)\n"
     if pp.homogeneous == "None"
         hom_params = Empty{FT}()
     elseif pp.homogeneous == "ABHOM"
@@ -248,7 +247,7 @@ function run_parcel(IC, t_0, t_end, pp)
         throw("Unrecognized homogeneous mode")
     end
 
-    println("Condensation growth: ", pp.condensation_growth)
+    info *= "Condensation growth: $(pp.condensation_growth)\n"
     if pp.condensation_growth == "None"
         ce_params = Empty{FT}()
     elseif pp.condensation_growth == "Condensation"
@@ -257,15 +256,15 @@ function run_parcel(IC, t_0, t_end, pp)
         throw("Unrecognized condensation growth mode")
     end
 
-    println("Deposition growth: ", pp.deposition_growth)
+    info *= "Deposition growth: $(pp.deposition_growth)\n"
     if pp.deposition_growth == "None"
         ds_params = Empty{FT}()
     elseif pp.deposition_growth == "Deposition"
         ds_params = DepParams{FT}(pp.aps, pp.tps)
     else
         throw("Unrecognized deposition growth mode")
     end
-    println(" ")
+    @info info
 
     # Parameters for the ODE solver
     p = (

parcel/ParcelTendencies.jl

@@ -141,9 +141,7 @@ function INPC_model(params::Frostenberg_stochastic, state)
     μ = CMI_het.INP_concentration_mean(T)
     g = ip.σ * sqrt(2 * γ)
 
-    dln_INPC =
-        -γ * (ln_INPC - μ) * const_dt +
-        g * sqrt(const_dt) * rand(DS.Normal(0, 1))
+    dln_INPC = -γ * (ln_INPC - μ) * const_dt + g * sqrt(const_dt) * randn()
 
     return dln_INPC / const_dt
 end