diff --git a/.github/workflows/DocCleanUp.yml b/.github/workflows/DocCleanUp.yml new file mode 100644 index 00000000..f3d42b6b --- /dev/null +++ b/.github/workflows/DocCleanUp.yml @@ -0,0 +1,28 @@ +name: Doc Preview Cleanup + +on: + pull_request: + types: [closed] + +jobs: + doc-preview-cleanup: + runs-on: ubuntu-latest + steps: + - name: Checkout gh-pages branch + uses: actions/checkout@v2 + with: + ref: gh-pages + + - name: Delete preview and history + run: | + git config user.name "Documenter.jl" + git config user.email "documenter@juliadocs.github.io" + git rm -rf "previews/PR$PRNUM" + git commit -m "delete preview" + git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree}) + env: + PRNUM: ${{ github.event.number }} + + - name: Push changes + run: | + git push --force origin gh-pages-new:gh-pages \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index c3141064..d689095e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,5 +5,6 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" diff --git a/docs/ThreeDimensionalInput.jl b/docs/ThreeDimensionalInput.jl deleted file mode 100644 index 0e77944a..00000000 --- a/docs/ThreeDimensionalInput.jl +++ /dev/null @@ -1,95 +0,0 @@ -using Thermodynamics -using Thermodynamics.TemperatureProfiles -using Thermodynamics.TestedProfiles -using UnPack -using CLIMAParameters -using CLIMAParameters.Planet -using Plots -import Thermodynamics -Thermodynamics.print_warning() = false - -struct EarthParameterSet <: AbstractEarthParameterSet end; -const param_set = EarthParameterSet(); -FT = Float64; -thermo_dir = dirname(dirname(pathof(Thermodynamics))); -include(joinpath(thermo_dir, "docs", "plot_helpers.jl")); -profiles = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT}); -@unpack ρ, e_int, q_tot = profiles - -dims = (10, 10, 10); -ρ = range(min(ρ...), stop = max(ρ...), length = dims[1]); -e_int = range(min(e_int...), stop = max(e_int...), length = dims[2]); -q_tot = range(min(q_tot...), stop = max(q_tot...), length = dims[3]); - -ρ_all = Array{FT}(undef, prod(dims)); -e_int_all = Array{FT}(undef, prod(dims)); -q_tot_all = Array{FT}(undef, prod(dims)); - -linear_indices = LinearIndices((1:dims[1], 1:dims[2], 1:dims[3])); -TS = Array{Union{ThermodynamicState, Nothing}}(undef, prod(dims)); -TS_no_err = Array{ThermodynamicState}(undef, prod(dims)); - -@inbounds for i in linear_indices.indices[1] - @inbounds for j in linear_indices.indices[2] - @inbounds for k in linear_indices.indices[3] - p = linear_indices[i, j, k] - ρ_all[p] = ρ[i] - e_int_all[p] = e_int[j] - q_tot_all[p] = q_tot[k] - - Thermodynamics.error_on_non_convergence() = false - TS_no_err[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k]) - Thermodynamics.error_on_non_convergence() = true - # @show p/prod(linear_indices.indices)*100 - try - TS[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k]) - catch - TS[p] = nothing - end - end - end -end - -# Full 3D scatter plot -function save_masked_plot3D(TS_no_err, mask, title, filename) - ρ_mask = ρ_all[mask] - T_mask = air_temperature.(TS_no_err[mask]) - q_tot_mask = q_tot_all[mask] - pts = (ρ_mask, T_mask, q_tot_mask) - Plots.plot(pts..., seriestype = :scatter, markersize = 7) - plot!( - xlabel = "Density", - ylabel = "Temperature", - zlabel = "Total specific humidity", - title = "$title", - ) - savefig(filename) -end; - -save_masked_plot3D(TS_no_err, TS .== nothing, "NC", "Scatter3DNonConverged.svg"); -save_masked_plot3D(TS_no_err, TS .!= nothing, "C", "Scatter3DConverged.svg"); - -# 2D binned scatter plots -function save_masked_plot2D_slices(TS_no_err, mask, title, filename) - ρ_mask = ρ_all[mask] - T_mask = air_temperature.(TS_no_err[mask]) - q_tot_mask = q_tot_all[mask] - save_binned_surface_plots(ρ_mask, T_mask, q_tot_mask, title, filename) -end; - -save_masked_plot2D_slices( - TS_no_err, - TS .== nothing, - "NC", - "Slices2DNonConverged.svg", -); -save_masked_plot2D_slices( - TS_no_err, - TS .!= nothing, - "C", - "Slices2DConverged.svg", -); - -@warn "Note that the temperature axis for the non-converged -plot is not necessarily accurate, since the temperatures are -the result of a non-converged saturation adjustment" diff --git a/docs/src/DevDocs.md b/docs/src/DevDocs.md index acd1db9d..038b3474 100644 --- a/docs/src/DevDocs.md +++ b/docs/src/DevDocs.md @@ -1,4 +1,7 @@ -# Input space exploration +# Saturation adjustment input space convergence maps + +The saturation adjustment procedure requires solving a non-linear +equation. In the [Tested Profiles](@ref) section, we plotted the tested thermodynamic states. In this section, we explore the convergence of the @@ -11,18 +14,35 @@ or likely to be observed in climate simulations, but showing the convergence space helps illustrate the buffer between our tested profiles and the nearest space where convergence fails. +This section is dedicated to monitoring the status and improvement +of the performance and robustness of various numerical methods +in solving the saturation adjustment equations for various thermodynamic +formulations. + +!!! note + + `dims` in `docs/src/saturation_adjustment.jl` is currently set to + ``dims = (6, 6, 6);`` to avoid heavy computations in the doc build, but you + may want to increase it to, e.g., ``dims = (10, 10, 10);`` when running locally + to see a higher resolution map. + ```@example -include(joinpath(@__DIR__, "..", "ThreeDimensionalInput.jl")) +include("saturation_adjustment.jl") ``` -## Converged cases (3D view) -![](Scatter3DConverged.svg) - -## Non-converged cases (3D view) -![](Scatter3DNonConverged.svg) +## 3D space +| Numerical method | Converged | Non-converged | +:-----------------:|:-----------------:|:---------------------: +SecantMethod | ![](3DSpace_converged_SecantMethod.svg) | ![](3DSpace_non_converged_SecantMethod.svg) +NewtonsMethod | ![](3DSpace_converged_NewtonsMethod.svg) | ![](3DSpace_non_converged_NewtonsMethod.svg) +NewtonsMethodAD | ![](3DSpace_converged_NewtonsMethodAD.svg) | ![](3DSpace_non_converged_NewtonsMethodAD.svg) +RegulaFalsiMethod | ![](3DSpace_converged_RegulaFalsiMethod.svg) | ![](3DSpace_non_converged_RegulaFalsiMethod.svg) -## Converged cases (2D view), binned by total specific humidity -![](Slices2DConverged.svg) +## 2D slices, binned by total specific humidity -## Non-converged cases (2D view), binned by total specific humidity -![](Slices2DNonConverged.svg) +| Numerical method | Converged | Non-converged | +:-----------------:|:-----------------:|:---------------------: +SecantMethod | ![](2DSlice_converged_SecantMethod.svg) | ![](2DSlice_non_converged_SecantMethod.svg) +NewtonsMethod | ![](2DSlice_converged_NewtonsMethod.svg) | ![](2DSlice_non_converged_NewtonsMethod.svg) +NewtonsMethodAD | ![](2DSlice_converged_NewtonsMethodAD.svg) | ![](2DSlice_non_converged_NewtonsMethodAD.svg) +RegulaFalsiMethod | ![](2DSlice_converged_RegulaFalsiMethod.svg) | ![](2DSlice_non_converged_RegulaFalsiMethod.svg) diff --git a/docs/src/saturation_adjustment.jl b/docs/src/saturation_adjustment.jl new file mode 100644 index 00000000..2b92aa6a --- /dev/null +++ b/docs/src/saturation_adjustment.jl @@ -0,0 +1,242 @@ +import Thermodynamics +import Thermodynamics.TemperatureProfiles + +const TD = Thermodynamics +const TP = TD.TemperatureProfiles + +import RootSolvers +const RS = RootSolvers + +import CLIMAParameters +const CP = CLIMAParameters + +import UnPack +import Plots +TD.print_warning() = false + +struct EarthParameterSet <: CP.AbstractEarthParameterSet end; +const param_set = EarthParameterSet(); +FT = Float64; + +const src_dir = dirname(dirname(pathof(Thermodynamics))); +profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array{FT}); +UnPack.@unpack ρ, q_tot = profiles +T_true = profiles.T +prof_pts = (ρ, T_true, q_tot) + +dims = (6, 6, 6); +ρ = range(min(ρ...), stop = max(ρ...), length = dims[1]); +T_true = range(min(T_true...), stop = max(T_true...), length = dims[2]); +q_tot = range(min(q_tot...), stop = max(q_tot...), length = dims[3]); + +ρ_all = Array{FT}(undef, prod(dims)); +T_true_all = Array{FT}(undef, prod(dims)); +q_tot_all = Array{FT}(undef, prod(dims)); + +linear_indices = LinearIndices((1:dims[1], 1:dims[2], 1:dims[3])); + +numerical_methods = ( + RS.SecantMethod, + RS.NewtonsMethod, + # RS.NewtonsMethodAD, # we need to relax diagonalization somewhere + # RS.RegulaFalsiMethod, # hit assertion error, bounds need adjusted for 3D space +) + +ts = Dict( + NM => Array{Union{TD.ThermodynamicState, Nothing}}(undef, prod(dims)) + for NM in numerical_methods +) +ts_no_err = Dict( + NM => Array{TD.ThermodynamicState}(undef, prod(dims)) + for NM in numerical_methods +) + +@inbounds for i in linear_indices.indices[1] + @inbounds for j in linear_indices.indices[2] + @inbounds for k in linear_indices.indices[3] + n = linear_indices[i, j, k] + ρ_all[n] = ρ[i] + T_true_all[n] = T_true[j] + q_tot_all[n] = q_tot[k] + e_int = TD.internal_energy( + param_set, + T_true[j], + TD.PhasePartition(q_tot[k]), + ) + + @inbounds for NM in numerical_methods + TD.error_on_non_convergence() = false + ts_no_err[NM][n] = TD.PhaseEquil_dev_only( + param_set, + ρ[i], + e_int, + q_tot[k]; + sat_adjust_method = NM, + maxiter = 10, + ) + TD.error_on_non_convergence() = true + # @show n / prod(dims) * 100 + try + ts[NM][n] = TD.PhaseEquil_dev_only( + param_set, + ρ[i], + e_int, + q_tot[k]; + sat_adjust_method = NM, + maxiter = 10, + ) + catch + ts[NM][n] = nothing + end + end + end + end +end + +# folder = "sat_adjust_analysis" +folder = @__DIR__ +mkpath(folder) + +function save_binned_surface_plots( + x, + y, + z, + title, + filename, + n_plots = (3, 3), + z_label_prefix = "z", + n_digits = 5; + xlims = (:auto, :auto), + ylims = (:auto, :auto), + label = label, + ref_points = nothing, +) + n_z_partitions = prod(n_plots) + local z_min_global, z_max_global + try + z_min_global = min(z...) + z_max_global = max(z...) + catch + z_min_global = 0 + z_max_global = 1 + end + Δz = (z_max_global - z_min_global) / n_z_partitions + z_min = ntuple(i -> z_min_global + (i - 1) * Δz, n_z_partitions) + z_max = ntuple(i -> z_min_global + (i) * Δz, n_z_partitions) + p = [] + for i in 1:n_z_partitions + mask = z_min[i] .<= z .<= z_max[i] + x_i = x[mask] + y_i = y[mask] + sz_min = string(z_min[i])[1:min(n_digits, length(string(z_min[i])))] + sz_max = string(z_max[i])[1:min(n_digits, length(string(z_max[i])))] + p_i = Plots.plot( + x_i, + y_i, + title = "$(title), in ($sz_min, $sz_max)", + seriestype = :scatter, + markersize = 5, + label = label, + xlims = xlims, + ylims = ylims, + ) + if ref_points ≠ nothing + ref_mask = z_min[i] .<= ref_points[3] .<= z_max[i] + Plots.plot!( + ref_points[1][ref_mask], + ref_points[2][ref_mask], + seriestype = :scatter, + markersize = 5, + label = "ref points", + ) + end + push!(p, p_i) + end + Plots.plot(p..., layout = n_plots, legend = false) + Plots.savefig(filename) +end; +# Full 3D scatter plot +function plot3D(ts_no_err, ts, NM; converged) + mask = converged ? ts .≠ nothing : ts .== nothing + + c_name = converged ? "converged" : "non_converged" + label = converged ? "converged" : "non-converged" + casename = converged ? "converged" : "non-converged" + nm_name = nameof(NM) + filename = "3DSpace_$(c_name)_$nm_name.svg" + + ρ_mask = ρ_all[mask] + q_tot_mask = q_tot_all[mask] + T_mask = T_true_all[mask] + pts = (ρ_mask, T_mask, q_tot_mask) + Plots.plot( + pts..., + color = "blue", + seriestype = :scatter, + markersize = 7, + label = casename, + ) + Plots.plot!( + prof_pts..., + color = "red", + seriestype = :scatter, + markersize = 7, + label = "tested thermo profiles", + ) + Plots.plot!( + xlabel = "Density", + ylabel = "Temperature", + zlabel = "Total specific humidity", + title = "3D input to PhaseEquil", + xlims = (min(ρ_all...), max(ρ_all...)), + ylims = (min(T_true_all...), max(T_true_all...)), + zlims = (min(q_tot_all...), max(q_tot_all...)), + camera = (50, 50), + # camera = (50,70), + ) + Plots.savefig(joinpath(folder, filename)) +end + +# 2D binned scatter plots +function plot2D_slices(ts_no_err, ts, NM; converged) + mask = converged ? ts .≠ nothing : ts .== nothing + ρ_mask = ρ_all[mask] + T_mask = T_true_all[mask] + q_tot_mask = q_tot_all[mask] + c_name = converged ? "converged" : "non_converged" + label = converged ? "converged" : "non-converged" + short_name = converged ? "C" : "NC" + nm_name = nameof(NM) + filename = "2DSlice_$(c_name)_$nm_name.svg" + filename = joinpath(folder, filename) + save_binned_surface_plots( + ρ_mask, + T_mask, + q_tot_mask, + short_name, + filename; + xlims = (min(ρ_all...), max(ρ_all...)), + ylims = (min(T_true_all...), max(T_true_all...)), + label = label, + ref_points = prof_pts, + ) +end + +# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846 +ENV["GKSwstype"] = "nul" + +for NM in numerical_methods + plot3D(ts_no_err[NM], ts[NM], NM; converged = false) + plot3D(ts_no_err[NM], ts[NM], NM; converged = true) + plot2D_slices(ts_no_err[NM], ts[NM], NM; converged = true) + plot2D_slices(ts_no_err[NM], ts[NM], NM; converged = false) +end + +convergence_percent = Dict() +for NM in numerical_methods + convergence_percent[NM] = count(ts[NM] .≠ nothing) / length(ts[NM]) +end +println("Convergence percentages:") +for (k, v) in convergence_percent + println("$k = $v") +end