diff --git a/src/RayleighBenardConvection/Convection.dat b/src/RayleighBenardConvection/Convection.dat
index 7a40b85..6b56f0e 100644
--- a/src/RayleighBenardConvection/Convection.dat
+++ b/src/RayleighBenardConvection/Convection.dat
@@ -124,6 +124,7 @@
out_j2_dev_stress = 1
out_j2_strain_rate = 1
out_temperature = 1
+ out_visc_total = 1
# AVD phase viewer output options (requires activation)
diff --git a/src/RayleighBenardConvection/Convection_Dash.jl b/src/RayleighBenardConvection/Convection_Dash.jl
index cc23fb7..8fc0b4d 100644
--- a/src/RayleighBenardConvection/Convection_Dash.jl
+++ b/src/RayleighBenardConvection/Convection_Dash.jl
@@ -44,14 +44,17 @@ function convection(; host=HTTP.Sockets.localhost, port=8050, width="80vw", heig
make_title(title_app),
dbc_row([
dbc_col([
- make_plot("",cmaps, width=width, height=height), # show graph
- make_plot_controls(), # show media buttons
- make_id_label(), # show user id
+ make_plot("",cmaps, width=width, height=height),
+ dbc_collapse(make_nu_ra_plot(width=width), id="collapse-nu-ra", is_open=true),
+ dbc_collapse(make_temp_profile_plot(width=width), id="collapse-temp-profile", is_open=false),
+ make_plot_controls(),
+ make_id_label(),
]),
dbc_col([
- make_time_card(), # show simulation time info
- make_menu(cmaps, show_field="temperature [°C]"), # show menu with simulation parameters, rheological parameters, and plotting parameters
- make_run_button() # show the run simulation button
+ make_time_card(),
+ make_extra_plot_controls(),
+ make_menu(cmaps, show_field="temperature [°C]"),
+ make_run_button()
])
]),
@@ -63,6 +66,9 @@ function convection(; host=HTTP.Sockets.localhost, port=8050, width="80vw", heig
dcc_store(id="last_timestep", data="0"),
dcc_store(id="update_fig", data="0"),
+ # Accumulates Nu values and Ra for the Nu/Ra chart
+ dcc_store(id="nu_ra_store", data=Dict("timesteps"=>[], "times"=>[], "nu"=>[], "ra"=>0.0)),
+
# Start an interval that updates the number every second
dcc_interval(id="session-interval", interval=100, n_intervals=0, disabled=true)
])
@@ -259,7 +265,7 @@ function convection(; host=HTTP.Sockets.localhost, port=8050, width="80vw", heig
if has_pvd_file(OutFile, user_dir)
Timestep, _, Time = read_LaMEM_simulation(OutFile, user_dir) # all timesteps
id = findall(Timestep .== cur_t)[1]
- if trigger == "button-start.n_clicks" || trigger == "button-play.n_clicks"
+ if trigger == "button-start.n_clicks"
cur_t = 0
id = 1
elseif trigger == "button-last.n_clicks"
@@ -318,7 +324,129 @@ function convection(; host=HTTP.Sockets.localhost, port=8050, width="80vw", heig
return label_timestep, label_time, current_timestep, fig_cross, add_units(fields_available), add_units(fields_available)
end
- #
+ # Accumulate Nu values for each new timestep; derive Ra from simulation output fields
+ callback!(app,
+ Dash.Output("nu_ra_store", "data"),
+ Input("last_timestep", "data"),
+ Input("button-run", "n_clicks"),
+ State("nu_ra_store", "data"),
+ State("session-id", "data"),
+ prevent_initial_call=true
+ ) do last_timestep, n_run, store_data, session_id
+
+ trigger = get_trigger()
+
+ if trigger == "button-run.n_clicks"
+ return Dict("timesteps"=>[], "times"=>[], "nu"=>[], "ra"=>0.0)
+ end
+
+ ts_list = collect(Any, store_data["timesteps"])
+ time_list = collect(Any, store_data["times"])
+ nu_list = collect(Any, store_data["nu"])
+ ra_stored = store_data["ra"]
+ last_computed = isempty(ts_list) ? -1 : Int(ts_list[end])
+ last_t = parse(Int, last_timestep)
+
+ last_t <= last_computed && return store_data
+
+ user_dir = simulation_directory(session_id, clean=false)
+ !has_pvd_file(OutFile, user_dir) && return store_data
+
+ Timestep, _, Time = read_LaMEM_simulation(OutFile, user_dir)
+ Ra = ra_stored # Ra is fixed once latched from the first timestep
+
+ for i in eachindex(Timestep)
+ t = Timestep[i]
+ t > last_computed && t <= last_t || continue
+
+ # Nu needs only temperature — always attempt this
+ nu = try
+ raw, _ = read_LaMEM_timestep(OutFile, t, user_dir)
+
+ T = raw.fields.temperature[:, 1, :]'
+ z = raw.z.val[1, 1, :]
+ nu_val = compute_nusselt(T, z)
+
+ # Ra: needs density + visc_total; latch once, fail gracefully
+ if Ra == 0.0
+ try
+ rho = raw.fields.density[:, 1, :]'
+ visc = raw.fields.visc_total[:, 1, :]'
+ nx = size(T, 2)
+ ΔT_ = sum(T[1, :]) / nx - sum(T[end, :]) / nx
+ d_ = (z[end] - z[1]) * 1e3
+ ρ_ = sum(rho) / length(rho)
+ η_ = 10^(sum(visc) / length(visc))
+ if ΔT_ > 0
+ Ra = ρ_ * 9.81 * 3e-5 * ΔT_ * d_^3 / (η_ * 8.7e-7)
+ end
+ catch e_ra
+ println("Ra computation skipped at t=$t: ", e_ra)
+ end
+ end
+
+ nu_val
+ catch e_nu
+ println("Nu computation error at t=$t: ", e_nu)
+ NaN
+ end
+
+ if isfinite(nu)
+ push!(ts_list, t)
+ push!(time_list, Time[i])
+ push!(nu_list, nu)
+ end
+ end
+
+ return Dict("timesteps"=>ts_list, "times"=>time_list, "nu"=>nu_list, "ra"=>Ra)
+ end
+
+ # Render the Nu / Ra chart whenever the store is updated
+ callback!(app,
+ Dash.Output("figure_nu_ra", "figure"),
+ Input("nu_ra_store", "data"),
+ State("n_timesteps", "value"),
+ prevent_initial_call=true
+ ) do store_data, n_timesteps
+ return create_nu_ra_figure(store_data, n_timesteps)
+ end
+
+ # Toggle visibility of the extra diagnostic plots
+ callback!(app,
+ Dash.Output("collapse-nu-ra", "is_open"),
+ Dash.Output("collapse-temp-profile", "is_open"),
+ Input("switch-extra-plots", "value"),
+ prevent_initial_call=false
+ ) do switch_vals
+ show_nu = !isnothing(switch_vals) && "Nu/Ra over time" in switch_vals
+ show_T = !isnothing(switch_vals) && "Mean T(z) profile" in switch_vals
+ return show_nu, show_T
+ end
+
+ # Update the mean T(z) profile for the currently displayed timestep
+ callback!(app,
+ Dash.Output("figure_temp_profile", "figure"),
+ Input("current_timestep", "data"),
+ State("session-id", "data"),
+ State("switch-extra-plots", "value"),
+ prevent_initial_call=true
+ ) do current_timestep, session_id, switch_vals
+ empty_fig = (data=[], layout=(autosize=true, margin=attr(t=45, b=40, l=55, r=10)))
+ !isnothing(switch_vals) && "Mean T(z) profile" in switch_vals || return empty_fig
+
+ cur_t = parse(Int, current_timestep)
+ user_dir = simulation_directory(session_id, clean=false)
+ has_pvd_file(OutFile, user_dir) || return empty_fig
+
+ try
+ x, z, data2D, _, _ = get_data(OutFile, cur_t, "temperature [°C]", user_dir)
+ return create_temp_profile_figure(data2D, z)
+ catch
+ return empty_fig
+ end
+ end
+
+ #
callback!(app,
Dash.Output("contour_option", "disabled"),
Input("switch-contour", "value")) do switch_contour
diff --git a/src/RayleighBenardConvection/dash_functions_convection.jl b/src/RayleighBenardConvection/dash_functions_convection.jl
index 2bb94c2..8f333bb 100644
--- a/src/RayleighBenardConvection/dash_functions_convection.jl
+++ b/src/RayleighBenardConvection/dash_functions_convection.jl
@@ -2,6 +2,229 @@ using DelimitedFiles
make_geometry_parameters() = nothing
+"""
+ Nu = compute_nusselt(data2D, z)
+
+Compute the Nusselt number from a 2D temperature field.
+
+`data2D[iz, ix]` is temperature in °C, `iz=1` is the bottom boundary and
+`iz=end` is the top boundary. `z` is the 1-D coordinate vector (km),
+running from −H (bottom) to 0 (top).
+
+Nu = 1 for pure conduction; Nu > 1 indicates convective heat transport.
+"""
+function compute_nusselt(data2D, z)
+ nz, nx = size(data2D)
+
+ T_top = data2D[end, :] # top boundary (z = 0)
+ T_bot = data2D[1, :] # bottom boundary (z = -H)
+
+ mean_T_top = sum(T_top) / nx
+ mean_T_bot = sum(T_bot) / nx
+ ΔT_actual = mean_T_top - mean_T_bot # negative (hot bottom, cold top)
+
+ abs(ΔT_actual) < 1e-10 && return NaN
+
+ # dT/dz at the top boundary [°C / km]
+ Δz = z[end] - z[end-1] # > 0
+ dT_dz_top = (data2D[end, :] .- data2D[end-1, :]) ./ Δz # < 0
+ mean_dT_dz = sum(dT_dz_top) / nx
+
+ H = z[end] - z[1] # domain height, positive [km]
+
+ # Nu = H * mean(∂T/∂z|_top) / (T_top - T_bot)
+ # Both numerator and denominator are negative → Nu > 0.
+ # For pure conduction: ∂T/∂z = -ΔT/H → Nu = 1.
+ return H * mean_dT_dz / ΔT_actual
+end
+
+"""
+ nu, Ra = compute_nu_ra(data; α=3e-5, κ=8.7e-7)
+
+Compute the Nusselt and Rayleigh numbers from a single LaMEM `CartData`
+object (as returned by `read_LaMEM_timestep`).
+
+ρ and η are read directly from the `density` and `visc_total` output fields;
+ΔT and d are derived from the temperature field and domain coordinates.
+Only the thermal expansion coefficient α [1/K] and the thermal diffusivity
+κ [m²/s] must be assumed (exercise-sheet defaults: α = 3×10⁻⁵ K⁻¹,
+κ = k/(ρ cₚ) = 3/(3300·1050) ≈ 8.7×10⁻⁷ m²/s).
+
+g = 9.81 m/s² is read from the simulation setup.
+"""
+function compute_nu_ra(data; α=3e-5, κ=8.7e-7, g=9.81)
+ T = data.fields.temperature[:, 1, :]' # [nz, nx], °C
+ z = data.z.val[1, 1, :] # z coordinates, km (bottom → top)
+
+ nu = compute_nusselt(T, z)
+
+ # ΔT from actual boundary temperatures
+ nx = size(T, 2)
+ ΔT = sum(T[1, :]) / nx - sum(T[end, :]) / nx # T_bot − T_top [K]
+ ΔT <= 0 && return nu, 0.0
+
+ # Domain height (km → m)
+ d = (z[end] - z[1]) * 1e3
+
+ # Representative density: spatial mean [kg/m³]
+ rho = data.fields.density[:, 1, :]'
+ ρ = sum(rho) / length(rho)
+
+ # Representative viscosity: geometric mean via log-space average [Pa·s]
+ # visc_total is stored as log₁₀(Pa·s) in LaMEM output
+ visc = data.fields.visc_total[:, 1, :]'
+ η = 10^(sum(visc) / length(visc))
+
+ Ra = ρ * g * α * ΔT * d^3 / (η * κ)
+ return nu, Ra
+end
+
+"""
+ fig = create_nu_ra_figure(store_data, n_timesteps)
+
+Build a Plotly figure that shows the Nusselt number vs time and annotates
+the Rayleigh number. `store_data` is the Dict stored in `nu_ra_store`;
+`n_timesteps` is used to fix the x-axis upper bound.
+"""
+function create_nu_ra_figure(store_data, n_timesteps)
+ ts = store_data["timesteps"]
+ times = store_data["times"]
+ nu_vals = store_data["nu"]
+ ra = store_data["ra"]
+
+ ra_str = ra > 0 ? "Ra = $(round(ra, sigdigits=3))" : "Ra = --"
+
+ base_layout = (
+ title = attr(text="Nu vs Time | $ra_str", font=attr(size=13)),
+ xaxis = attr(title="Time [Myrs]", range=[0, nothing]),
+ yaxis = attr(title="Nu", rangemode="tozero"),
+ margin = attr(t=45, b=40, l=55, r=10),
+ autosize = true,
+ legend = attr(x=0.01, y=0.99, xanchor="left", yanchor="top",
+ bgcolor="rgba(255,255,255,0.6)"),
+ showlegend = true,
+ )
+
+ isempty(nu_vals) && return (data=[], layout=base_layout)
+
+ trace_nu = scatter(
+ x = times,
+ y = nu_vals,
+ mode = "lines+markers",
+ name = "Nu",
+ marker = attr(size=4),
+ line = attr(color="steelblue", width=2),
+ customdata = ts,
+ hovertemplate = "Time: %{x:.3f} Myrs
Timestep: %{customdata}
Nu: %{y:.3f}",
+ )
+
+ # reference line Nu = 1 (pure conduction)
+ trace_cond = scatter(
+ x = [times[1], times[end]],
+ y = [1.0, 1.0],
+ mode = "lines",
+ name = "Nu = 1 (conductive)",
+ line = attr(color="firebrick", dash="dash", width=1),
+ hoverinfo = "skip",
+ )
+
+ return (
+ data = [trace_nu, trace_cond],
+ layout = base_layout,
+ config = (edits = (shapePosition=true,),),
+ )
+end
+
+"""
+ make_nu_ra_plot(; width, height)
+
+Return a `dbc_row` containing the Nu-vs-time graph element.
+"""
+function make_nu_ra_plot(; width="80vw", height="22vh")
+ return dbc_row([
+ dcc_graph(
+ id = "figure_nu_ra",
+ figure = (data=[], layout=(title="Nu vs Time | Ra = --",
+ autosize=true,
+ margin=attr(t=45, b=40, l=55, r=10))),
+ style = attr(width=width, height=height),
+ )
+ ])
+end
+
+
+"""
+ create_temp_profile_figure(data2D, z)
+
+Build a Plotly figure of the horizontally averaged temperature profile.
+`data2D[iz, ix]` is temperature in °C; `z` is the depth coordinate vector [km].
+"""
+function create_temp_profile_figure(data2D, z)
+ nx = size(data2D, 2)
+ T_mean = vec(sum(data2D, dims=2)) ./ nx
+
+ trace = scatter(
+ x = T_mean,
+ y = z,
+ mode = "lines",
+ line = attr(color="steelblue", width=2),
+ hovertemplate = "T: %{x:.1f} °C
z: %{y:.1f} km",
+ showlegend = false,
+ )
+
+ return (
+ data = [trace],
+ layout = (
+ title = attr(text="Mean Temperature Profile", font=attr(size=13)),
+ xaxis = attr(title="Temperature [°C]"),
+ yaxis = attr(title="Depth [km]"),
+ margin = attr(t=45, b=40, l=55, r=10),
+ autosize = true,
+ ),
+ config = (edits = (shapePosition=true,),),
+ )
+end
+
+"""
+ make_temp_profile_plot(; width, height)
+
+Return a `dbc_row` containing the mean-T-profile graph element.
+"""
+function make_temp_profile_plot(; width="80vw", height="30vh")
+ return dbc_row([
+ dcc_graph(
+ id = "figure_temp_profile",
+ figure = (data=[], layout=(title="Mean Temperature Profile",
+ autosize=true,
+ margin=attr(t=45, b=40, l=55, r=10))),
+ style = attr(width=width, height=height),
+ )
+ ])
+end
+
+"""
+ make_extra_plot_controls()
+
+Card with toggle switches for the diagnostic plots below the main figure.
+"""
+function make_extra_plot_controls()
+ return dbc_row([
+ html_p(),
+ dbc_card([
+ dbc_label("Diagnostic plots:", size="sm"),
+ dbc_checklist(
+ options = ["Nu/Ra over time", "Mean T(z) profile"],
+ value = ["Nu/Ra over time"],
+ id = "switch-extra-plots",
+ switch = true,
+ ),
+ ],
+ color = "secondary",
+ class_name = "mx-auto col-11",
+ outline = true),
+ html_p(),
+ ])
+end
"""
Returns an accordion menu containing the simulation parameters.