BOMEX plots & diagnostics

[mmr0]

work in progress reproducing more of the figures (currently 4 & 5) from Siebesma et al 2023. Not sure if these will live in the examples, but I thought they might be useful for JOSS paper?

[glwagner]

I think its great to put it in the examples!

[bischtob]

Plots! Yes, having plots ready will accelerate the JOSS paper!

[glwagner]

excited to see these!

[giordano]

Suggested change

	tke = @at (Center, Center, Center) (u′² + v′² + w′²) / 2
	tke = @at (Center, Center, Center) (u′² + v′² + w′²) / 2

[giordano]

This comment should eventually be removed

[mmr0]

I haven't been able to get this to work. cloud_mask = qˡ .> qˡ_thresh gives
ERROR: LoadError: Scalar indexing is disallowed.

But I thought broadcasting '.>' was supposed to avoid scalar indexing?

[mmr0]

@glwagner any suggestions?

[mmr0]

maybe cloud_mask = findall(qˡ .> qˡ_thresh)

[giordano]

What are you exactly trying to achieve here? qˡ .> qˡ_thresh would give you an array of zeros and ones, is that what you want? Even if you go beyond the scalar indexing issue, Maximum on the line below is undefined, was that supposed to be maximum?

[mmr0]

Yes I think that's what I want - a 3d array of cloud presence/absence, which I then reduce to a 2d (x,y) array which tells me if there is any cloud in each vertical column. I think maximum will achieve that (although not Maximum as you point out -- thanks!). Then I can compute total cloud cover as a fractional area.

From siebesma 2003:

The total cloud cover is defined as the fraction of vertical columns that contain cloud water and is therefore identical to the cloud cover that would be observed ideally by satellite.

I'm not sure what their qˡ_thresh is though

[giordano]

What do you eventually want to do with cloud_fraction? Is it ok to convert it to a simple Array, or you want the corresponding field? At the moment this chunk of code is unused.

[mmr0]

I want to pass it to the output writer alongside the integrated TKE

[giordano]

Ok, based on https://clima.github.io/OceananigansDocumentation/stable/appendix/library/#Oceananigans.OutputWriters.JLD2Writer-Tuple%7BAny,%20Any%7D do I understand correctly that cloud_fraction needs to be a Field, and an Array wouldn't work?

[glwagner]

I think that what you want is an abstract operation — like
op = a > b. This can be output directly, or used to form another abstract operation. If it is output, it will be wrapped in a Field and recomputed each time the output is fetched.

[glwagner]

If you invoke broadcasting, you will make a one-time computation but it will not recompute every time the computation is requested. This is useful in principle for plotting or whatever, but not with an output writer. Also a separate issue is that broadcasting is a bit borked. To do it right, you may need to preform the result. But note, for post processing you could also use abstract operations. So generally it’s preferred to use an operation. This also means that if desired, the post processing code can be computed online with little additional code changes.

[mmr0]

Oh and here are the new figures:

[giordano]

@mmr0 sorry, I did some git surgery because this PR had loads of unrelated changes coming from other branches, which made review a lot harder. Locally you'll have to run the following commands:

git fetch --all --prune
git checkout mgr/bomex-turb-quantities
git reset --hard origin/mgr/bomex-turb-quantities

[giordano]

Ok, I just built the Literate example locally, the problem with the size is simply that there are multiple large images. In particular,

Breeze.jl/examples/bomex.jl

Line 474 in b5c1df2

fig

alone, added in this PR, is over 700 kB.

[giordano]

One workaround is to display the files written to disk, instead of the fig objects:

diff --git a/examples/bomex.jl b/examples/bomex.jl
index e0f710d..c4cea1c 100644
--- a/examples/bomex.jl
+++ b/examples/bomex.jl
@@ -400,9 +400,8 @@ text!(axuv, -8.5, 2200, text="solid: u\ndashed: v", fontsize=12)
 
 fig[0, :] = Label(fig, "BOMEX: Mean profile evolution (Siebesma et al., 2003)", fontsize=18, tellwidth=false)
 
-save("bomex_profiles.png", fig) #src
-fig
-
+save("bomex_profiles.png", fig)
+# ![](bomex_profiles.png)
 
 # 1 x 2 panel plot showing vertical velocity variance and tke
 w′²t = FieldTimeSeries(filename, "w′²")
@@ -433,8 +432,8 @@ xlims!(axw, 0, 0.3)
 
 fig[0, :] = Label(fig, "BOMEX: turbulent profile evolution (Siebesma et al., 2003)", fontsize=18, tellwidth=false)
 
-save("bomex_var_profiles.png", fig) #src
-fig
+save("bomex_var_profiles.png", fig)
+# ![](bomex_var_profiles.png)
 
 # 3 x 2 panel plot showing turbulent flux profiles
 w′qˡ′t = FieldTimeSeries(filename, "w′qˡ′")
@@ -470,8 +469,8 @@ axislegend(axwqt, position=:rt)
 
 fig[0, :] = Label(fig, "BOMEX: Turbulent flux profile evolution (Siebesma et al., 2003)", fontsize=18, tellwidth=false)
 
-save("bomex_turb_profiles.png", fig) #src
-fig
+save("bomex_turb_profiles.png", fig)
+# ![](bomex_turb_profiles.png)
 
 # The simulation shows the development of a cloudy boundary layer with:
 # - Warming of the subcloud layer from surface fluxes

This would make the generated HTML a lot smaller (1.7 MB -> 120 kB for me)

[giordano]

Suggested change

	default_colours = Makie.wong_colors();
	colors = [default_colours[mod1(i, length(default_colours))] for i in 1:Nt];
	default_colours = Makie.wong_colors()
	colors = [default_colours[mod1(i, length(default_colours))] for i in 1:Nt]

[mmr0]

@mmr0 sorry, I did some git surgery because this PR had loads of unrelated changes coming from other branches, which made review a lot harder. Locally you'll have to run the following commands:

git fetch --all --prune
git checkout mgr/bomex-turb-quantities
git reset --hard origin/mgr/bomex-turb-quantities

Copy - thanks! Sorry as you may gather I'm new to this so very much appreciate the clear instructions. Will have to do this tomorrow as my local branch sits on the supercomputer which is undergoing maintenance today 🥲