diff --git a/Project.toml b/Project.toml index bdb05b22e7..90b19a4131 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a" Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d" +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d" ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f" @@ -19,6 +20,7 @@ CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3" @@ -41,6 +43,7 @@ ArgParse = "1" ArtifactWrappers = "0.2" Artifacts = "1" AtmosphericProfilesLibrary = "0.1.7" +BlockArrays = "1" ClimaComms = "0.6.4" ClimaCore = "0.14.12" ClimaDiagnostics = "0.2.4" @@ -51,6 +54,7 @@ CloudMicrophysics = "0.22.3" Dates = "1" DiffEqBase = "6.145" FastGaussQuadrature = "0.5, 1" +ForwardDiff = "0.10.36" Insolation = "0.9.2" Interpolations = "0.15.1" LazyArtifacts = "1" diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 56abd972b9..a2531eb718 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -86,6 +86,15 @@ eisenstat_walker_forcing_alpha: jvp_step_adjustment: help: "Amount by which the step size of the forward difference approximation of the Jacobian-vector product in the Krylov method should be scaled (only used if `use_krylov_method` is `true`)" value: 1.0 +use_exact_jacobian: + help: "Whether to solve the linear system using the exact Jacobian, which is computed with forward-mode automatic differentiation (default is `false`)" + value: true # false +n_steps_update_exact_jacobian: + help: "Number of timesteps between updates to the exact Jacobian (default is 0, which corresponds to updating the exact Jacobian on each Newton iteration)" + value: 10 # 0 +debug_approximate_jacobian: + help: "Whether to check the accuracy of the approximate Jacobian against the exact Jacobian every `n_steps_update_exact_jacobian` timesteps (default is `false`)" + value: false # Radiation rad: help: "Radiation model [`nothing` (default), `gray`, `clearsky`, `allsky`, `allskywithclear`]" diff --git a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml index a7e29317ff..faea36866c 100644 --- a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml +++ b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml @@ -7,7 +7,7 @@ turbconv: "prognostic_edmfx" implicit_diffusion: true implicit_sgs_advection: true approximate_linear_solve_iters: 2 -max_newton_iters_ode: 3 +max_newton_iters_ode: 10 edmfx_upwinding: first_order edmfx_entr_model: "Generalized" edmfx_detr_model: "Generalized" diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 498ac7b324..35067d7da7 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -295,7 +295,7 @@ weakdeps = ["SparseArrays"] ChainRulesCoreSparseArraysExt = "SparseArrays" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.27.5" diff --git a/examples/Manifest.toml b/examples/Manifest.toml index d247f216bf..710f48519a 100644 --- a/examples/Manifest.toml +++ b/examples/Manifest.toml @@ -319,7 +319,7 @@ version = "0.5.7" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.27.5" diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 8777ee59ea..bb2971a2c5 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -145,10 +145,17 @@ if ClimaComms.iamroot(config.comms_ctx) include( joinpath(pkgdir(CA), "regression_tests", "self_reference_or_path.jl"), ) + + # TODO: Improve design of Diagnostics to make Jacobian plotting less clunky. + (; dict_writer) = simulation.output_writers + Yₜ_end = similar(integrator.u) + CA.implicit_tendency!(Yₜ_end, integrator.u, integrator.p, integrator.t) + @info "Plotting" path = self_reference_or_path() # __build__ path (not job path) if path == :self_reference - make_plots(Val(Symbol(reference_job_id)), simulation.output_dir) + paths = [simulation.output_dir] + make_plots(Val(Symbol(reference_job_id)), paths, dict_writer, Yₜ_end) else main_job_path = joinpath(path, reference_job_id) nc_dir = joinpath(main_job_path, "nc_files") @@ -170,11 +177,11 @@ if ClimaComms.iamroot(config.comms_ctx) end paths = if isempty(readdir(nc_dir)) - simulation.output_dir + [simulation.output_dir] else [simulation.output_dir, nc_dir] end - make_plots(Val(Symbol(reference_job_id)), paths) + make_plots(Val(Symbol(reference_job_id)), paths, dict_writer, Yₜ_end) end @info "Plotting done" diff --git a/perf/Manifest.toml b/perf/Manifest.toml index 052f4d32a1..f72b326092 100644 --- a/perf/Manifest.toml +++ b/perf/Manifest.toml @@ -330,7 +330,7 @@ version = "0.5.6" GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.27.5" diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index b058c2d71a..4256d32e8a 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -49,6 +49,7 @@ import ClimaCoreSpectra: power_spectrum_2d using Poppler_jll: pdfunite import Base.Filesystem +import LinearAlgebra: diag const days = 86400 @@ -71,19 +72,27 @@ function make_plots(sim, _) end """ - make_plots(sim, simulation_path::AbstractString) - make_plots(sim, simulation_paths::Iterable{AbstractString}) + make_plots(sim, simulation_paths, dict_writer, Yₜ_end) -Plot the corresponding sets for the given `sim`. +Plot the available data for the given `sim`. + +When `simulation_paths` contains one string, the data at that path is plotted. -When `simulation_path` is a string, use the data in the path. +When `simulations_paths` contains multiple strings, the data at those paths is +compared side-by-side. -When `simulation_paths` is a collection of strings, use all those paths and making -side-by-side comparisons. +All of the data stored in the `dict_writer` is also plotted by default. """ -function make_plots(sim, simulation_path::AbstractString) - paths = [simulation_path] - make_plots(sim, paths) +function make_plots(sim, simulation_paths, dict_writer, Yₜ_end) + for (time_series_name, time_series) in dict_writer.dict + filename = if startswith(time_series_name, "ejac") + "exact_jacobian" + else + error("Unrecognized time series name $time_series_name") + end + plot_jacobian(simulation_paths[1], time_series, Yₜ_end, filename) + end + make_plots(sim, simulation_paths) end # The contour plot functions in ClimaAnalysis work by finding the nearest slice available. @@ -469,6 +478,372 @@ function plot_spectrum_with_line!(grid_loc, spectrum; exponent = -3.0) return nothing end +function plot_jacobian(output_path, jacobian_time_series, Yₜ_end, filename) + # Take the transpose of each Jacobian matrix so that its rows are plotted + # along the y-axis and its columns are plotted along the x-axis. + time_series = + Dict(map(((t, matrix),) -> t => matrix', [jacobian_time_series...])) + + n_max_plot_times = 6 + n_available_times = length(time_series) + available_times = sort(collect(keys(time_series))) + plot_times = if n_available_times < n_max_plot_times + available_times + else + n_times_to_skip = fld(n_available_times, n_max_plot_times) + available_times[1:n_times_to_skip:(n_max_plot_times * n_times_to_skip)] + end + + field_names = collect(CA.scalar_field_names(Yₜ_end)) + first_tick_positions = map(field_names) do field_name + tick_position = 1 + for prev_field_name in field_names + field_name == prev_field_name && return tick_position + prev_field = CA.MatrixFields.get_field(Yₜ_end, prev_field_name) + if prev_field isa CA.Fields.SpectralElementField + tick_position += 1 + else + tick_position += CA.Spaces.nlevels(axes(prev_field)) + end + end + end + last_field = CA.MatrixFields.get_field(Yₜ_end, field_names[end]) + last_index = if last_field isa CA.Fields.SpectralElementField + first_tick_positions[end] + else + first_tick_positions[end] + CA.Spaces.nlevels(axes(last_field)) - 1 + end + last_tick_positions = [(first_tick_positions[2:end] .- 1)..., last_index] + + center_tick_positions = (first_tick_positions .+ last_tick_positions) ./ 2 + boundary_tick_positions = [0.5, (last_tick_positions .+ 0.5)...] + index_ranges = map(UnitRange, first_tick_positions, last_tick_positions) + + limit_padding = length(index_ranges[1]) / 20 + limit_positions = + extrema(boundary_tick_positions) .+ (-limit_padding, limit_padding) + + tick_labels = map(field_names) do field_name + replace( + string(field_name), + "@name(" => "", + ".components.data" => "", + ":(" => "", + ')' => "", + ) + end + axis_kwargs = (; + limits = (limit_positions, limit_positions), + titlegap = 100, + yreversed = true, + xlabel = "∂Y index", + ylabel = "∂Yₜ index", + xlabelpadding = 60, + ylabelpadding = 60, + xticks = (center_tick_positions, tick_labels), + yticks = (center_tick_positions, tick_labels), + xticksvisible = false, + yticksvisible = false, + xticklabelrotation = pi / 4, + xticklabelpad = 50, + yticklabelpad = 80, + xminorticks = boundary_tick_positions, + yminorticks = boundary_tick_positions, + xminorticksvisible = true, + yminorticksvisible = true, + xminorticksize = 50, + yminorticksize = 50, + xminortickwidth = 10, + yminortickwidth = 10, + xgridvisible = false, + ygridvisible = false, + xminorgridvisible = true, + yminorgridvisible = true, + xminorgridwidth = 10, + yminorgridwidth = 10, + spinewidth = 10, + ) # Flip the y-axis so that the diagonal runs from top-left to bottom-right. + colorbar_kwargs = (; + size = 150, + labelpadding = 60, + ticklabelpad = 50, + ticksize = 50, + tickwidth = 10, + spinewidth = 10, + ) + figure_kwargs = (; figure_padding = 300, fontsize = 150) + plot_size = 5000 + figure_size = (2.2, cld(length(plot_times), 2)) .* plot_size + + typical_tendency_values = Dict( + CA.MatrixFields.@name(c.ρ) => 1e-8, + CA.MatrixFields.@name(c.uₕ.components.data.:1) => 1e-6, + CA.MatrixFields.@name(c.uₕ.components.data.:2) => 1e-6, + CA.MatrixFields.@name(c.ρe_tot) => 1e-1, + CA.MatrixFields.@name(c.ρq_tot) => 1e-7, + CA.MatrixFields.@name(c.ρq_liq) => 1e-8, + CA.MatrixFields.@name(c.ρq_ice) => 1e-8, + CA.MatrixFields.@name(c.ρq_rai) => 1e-8, + CA.MatrixFields.@name(c.ρq_sno) => 1e-8, + CA.MatrixFields.@name(c.sgs⁰.ρatke) => 1e-3, + CA.MatrixFields.@name(c.sgsʲs.:(1).ρa) => 1e-5, + CA.MatrixFields.@name(c.sgsʲs.:(1).mse) => 1e-1, + CA.MatrixFields.@name(c.sgsʲs.:(1).q_tot) => 1e-7, + CA.MatrixFields.@name(f.u₃.components.data.:1) => 1e-2, + CA.MatrixFields.@name(f.sgsʲs.:(1).u₃.components.data.:1) => 1e-2, + CA.MatrixFields.@name(sfc.T) => 1e-3, + CA.MatrixFields.@name(sfc.water) => 1e-7, + ) + rescaling_vector = ClimaComms.allowscalar( + Vector, + ClimaComms.device(Yₜ_end), + CA.Fields.column(Yₜ_end, 1, 1, 1), + ) + for (field_name, index_range) in zip(field_names, index_ranges) + block_max = CA.smooth_maximum(rescaling_vector[index_range]) + rescaling_vector[index_range] .= + block_max < 1e-10 ? typical_tendency_values[field_name] : block_max + end + rescaling_matrix = (inv.(rescaling_vector) * rescaling_vector')' + # Take the transpose of the rescaling matrix so that its rows and columns + # are consistent with the Jacobian matrices. + + sign_or_nan(x) = iszero(x) ? typeof(x)(NaN) : sign(x) + logabs_or_nan(x) = iszero(x) ? typeof(x)(NaN) : log10(abs(x)) + entry_sign(matrix) = sign_or_nan.(matrix) + entry_logabs(matrix) = logabs_or_nan.(matrix) + entry_logabs_rescaled(matrix) = entry_logabs(rescaling_matrix .* matrix) + function block_logabs_rescaled(matrix) + block_max_matrix = similar(matrix) + for row_index_range in index_ranges, col_index_range in index_ranges + block_max_matrix[row_index_range, col_index_range] .= + CA.smooth_maximum(matrix[row_index_range, col_index_range]) + end + return entry_logabs_rescaled(block_max_matrix) + end + function row_logabs_rescaled(matrix) + row_max_matrix = similar(matrix) + for row_index_range in index_ranges, col in axes(matrix, 1) + row_max_matrix[row_index_range, col] .= + CA.smooth_maximum(matrix[row_index_range, col]) + end + return entry_logabs_rescaled(row_max_matrix) + end + function block_bandwidth(matrix) + block_bandwidth_matrix = similar(matrix, Int) + for row_index_range in index_ranges, col_index_range in index_ranges + block = matrix[row_index_range, col_index_range] + bands = (1 - length(row_index_range)):(length(col_index_range) - 1) + block_bandwidth_matrix[row_index_range, col_index_range] .= + all(iszero, block) ? 0 : + findlast(band -> any(!iszero, diag(block, band)), bands) - + findfirst(band -> any(!iszero, diag(block, band)), bands) + 1 + end + return block_bandwidth_matrix + end + + bandwidth_matrices = map(t -> block_bandwidth(time_series[t]), plot_times) + sign_matrices = map(t -> entry_sign(time_series[t]), plot_times) + max_bandwidth = maximum(maximum, bandwidth_matrices) + + categorical(colormap, n) = CairoMakie.cgrad(colormap, n; categorical = true) + main_colormap = categorical(:tol_iridescent, 21) + bandwidth_colors = [CairoMakie.RGB(1, 1, 1), main_colormap[1:(end - 1)]...] + bandwidth_colormap = categorical(bandwidth_colors, max_bandwidth + 1) + sign_colormap = categorical(:RdBu_5, 2) + rescaling_colors = + setindex!(CairoMakie.to_colormap(:RdBu_11), CairoMakie.RGB(1, 1, 1), 6) + rescaling_colormap = categorical(rescaling_colors, 21) + + first_column_field = CA.Fields.column(Yₜ_end.c, 1, 1, 1) + coord_field = + CA.Fields.coordinate_field(CA.Fields.level(first_column_field, 1)) + coord = + ClimaComms.allowscalar(getindex, ClimaComms.device(Yₜ_end), coord_field) + column_loc_str = if coord isa CA.Geometry.XZPoint + "x = $(coord.x) m" + elseif coord isa CA.Geometry.XYZPoint + "x = $(coord.x) m, y = $(coord.y) m" + elseif coord isa CA.Geometry.LatLongZPoint + "lat = $(coord.lat)°, long = $(coord.long)°" + else + error("Unrecognized coordinate type $(typeof(coord))") + end + + page_file_paths = map(n -> joinpath(output_path, "$filename$n.pdf"), 1:6) + full_file_path = joinpath(output_path, "$filename.pdf") + + figure = CairoMakie.Figure(; size = figure_size, figure_kwargs...) + for (index, (t, bandwidth_matrix, sign_matrix)) in + enumerate(zip(plot_times, bandwidth_matrices, sign_matrices)) + grid_position = figure[cld(index, 2), (index - 1) % 2 + 1] + grid_layout = CairoMakie.GridLayout(grid_position) + title = "∂Yₜ/∂Y at $column_loc_str, t = $(CA.time_and_units_str(t))" + axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...) + bandwidth_plot = CairoMakie.heatmap!( + axis, + boundary_tick_positions, + boundary_tick_positions, + bandwidth_matrix[first_tick_positions, first_tick_positions]; + colormap = bandwidth_colormap, + colorrange = (-0.5, max_bandwidth + 0.5), + ) + CairoMakie.translate!(bandwidth_plot, 0, 0, -100) # Put plot under grid. + sign_plot = CairoMakie.heatmap!( + axis, + sign_matrix; + colormap = sign_colormap, + colorrange = (-1, 1), + ) + grid_sublayout = CairoMakie.GridLayout(grid_layout[1, 2]) + CairoMakie.Colorbar( + grid_sublayout[1, 1], + bandwidth_plot; + label = "Bandwidth of matrix block", + colorbar_kwargs..., + ) + CairoMakie.Colorbar( + grid_sublayout[2, 1], + sign_plot; + label = "Sign of matrix entry", + ticks = ([-0.5, 0.5], Makie.latexstring.(["-", "+"])), + colorbar_kwargs..., + labelpadding = 30, + ) + CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05)) + CairoMakie.rowgap!(grid_sublayout, CairoMakie.Relative(0.05)) + end + CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1)) + CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1)) + CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.08)) + CairoMakie.rowgap!(figure.layout, CairoMakie.Relative(0.05)) + CairoMakie.save(page_file_paths[1], figure) + + tickformat = Base.Fix1(map, x -> Makie.latexstring("10^{$(round(Int, x))}")) + for (index, (logabs_transform, value_type_str, is_rescaled)) in enumerate(( + (block_logabs_rescaled, "block", true), + (row_logabs_rescaled, "block row", true), + (entry_logabs_rescaled, "entry", true), + (entry_logabs, "entry", false), + )) + matrices = map(t -> logabs_transform(time_series[t]), plot_times) + max_logabs = + maximum(matrix -> maximum(filter(!isnan, matrix)), matrices) + figure = CairoMakie.Figure(; size = figure_size, figure_kwargs...) + for (index, (t, matrix)) in enumerate(zip(plot_times, matrices)) + grid_position = figure[cld(index, 2), (index - 1) % 2 + 1] + grid_layout = CairoMakie.GridLayout(grid_position) + title = "∂Yₜ/∂Y at $column_loc_str, t = $(CA.time_and_units_str(t))" + axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...) + plot_args = if value_type_str == "block" + submatrix = matrix[first_tick_positions, first_tick_positions] + (boundary_tick_positions, boundary_tick_positions, submatrix) + elseif value_type_str == "block row" + submatrix = matrix[first_tick_positions, :] + (boundary_tick_positions, 1:last_index, submatrix) + else + (matrix,) + end + plot = CairoMakie.heatmap!( + axis, + plot_args...; + colormap = main_colormap, + colorrange = (max_logabs - 6, max_logabs), + lowclip = main_colormap[1], + ) # Only color in the top 6 orders of magnitude. + CairoMakie.translate!(plot, 0, 0, -100) # Put plot under grid. + is_max = value_type_str == "entry" + label = "Absolute value of $(is_max ? " max of" : "") matrix \ + $value_type_str$(is_rescaled ? " [s⁻¹]" : "")" + CairoMakie.Colorbar( + grid_layout[1, 2], + plot; + label, + tickformat, + colorbar_kwargs..., + ) + CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05)) + end + CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1)) + CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1)) + CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.08)) + CairoMakie.rowgap!(figure.layout, CairoMakie.Relative(0.05)) + CairoMakie.save(page_file_paths[index + 1], figure) + end + + figure = CairoMakie.Figure(; size = (2.2, 1) .* plot_size, figure_kwargs...) + + grid_position = figure[1, 1] + title = "Rescaling vector for Yₜ" + ylimits = (minimum(rescaling_vector) / 2, maximum(rescaling_vector) * 2) + vector_axis_kwargs = (; + limits = (limit_positions, ylimits), + titlegap = 100, + yscale = log10, + xlabel = "Yₜ index", + ylabel = "Vector entry", + xlabelpadding = 60, + ylabelpadding = 60, + xticks = (center_tick_positions, tick_labels), + xticksvisible = false, + xticklabelrotation = pi / 4, + xticklabelpad = 50, + yticklabelpad = 50, + xminorticks = boundary_tick_positions, + xminorticksvisible = true, + xminorticksize = 50, + yticksize = 50, + xminortickwidth = 10, + ytickwidth = 10, + xgridvisible = false, + xminorgridvisible = true, + xminorgridwidth = 10, + ygridwidth = 10, + spinewidth = 10, + ) + axis = CairoMakie.Axis(grid_position; title, vector_axis_kwargs...) + CairoMakie.lines!( + axis, + 1:last_index, + rescaling_vector; + color = rescaling_colormap[end], + linewidth = 40, + ) + + matrix = logabs_or_nan.(rescaling_matrix) + grid_layout = CairoMakie.GridLayout(figure[1, 2]) + title = "Rescaling matrix for ∂Yₜ/∂Y" + axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...) + plot = CairoMakie.heatmap!( + axis, + boundary_tick_positions, + boundary_tick_positions, + matrix[first_tick_positions, first_tick_positions]; + colormap = rescaling_colormap, + ) + CairoMakie.translate!(plot, 0, 0, -100) # Put plot under grid. + CairoMakie.Colorbar( + grid_layout[1, 2], + plot; + label = "Matrix entry", + tickformat, + colorbar_kwargs..., + ) + CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05)) + + CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1)) + CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1)) + CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.15)) + CairoMakie.save(page_file_paths[end], figure) + + pdfunite() do unite + run(Cmd([unite, page_file_paths..., full_file_path])) + end + for page_file_path in page_file_paths + Filesystem.rm(page_file_path; force = true) + end +end + ColumnPlots = Union{ Val{:single_column_hydrostatic_balance_ft64}, Val{:single_column_radiative_equilibrium_gray}, diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl index b763073bad..c01ebd1d98 100644 --- a/regression_tests/ref_counter.jl +++ b/regression_tests/ref_counter.jl @@ -1,4 +1,4 @@ -179 +180 #= 179: diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index d0c6cbda5b..1523c2aa9d 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -48,6 +48,10 @@ include(joinpath("prognostic_equations", "zero_tendency.jl")) include(joinpath("prognostic_equations", "implicit", "implicit_tendency.jl")) include(joinpath("prognostic_equations", "implicit", "implicit_solver.jl")) +include(joinpath("prognostic_equations", "implicit", "approx_jacobian.jl")) +include(joinpath("prognostic_equations", "implicit", "exact_jacobian.jl")) +include(joinpath("prognostic_equations", "implicit", "debug_jacobian.jl")) +include(joinpath("prognostic_equations", "implicit", "dual_fixes.jl")) include(joinpath("prognostic_equations", "remaining_tendency.jl")) include(joinpath("prognostic_equations", "forcing", "large_scale_advection.jl")) # TODO: should this be in tendencies/? diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index 77ee03de1b..999dc2eea0 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -564,7 +564,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ρ_prev_level ) / local_geometry_level.J / ρ_level - @. detrʲ_prev_level = detrainment( + @. detrʲ_prev_level = detrainment_from_thermo_state( params, z_prev_level, z_sfc_halflevel, diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index c109e80df4..1903c5c5aa 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -4,6 +4,43 @@ import Thermodynamics as TD import ClimaCore: Spaces, Fields +function implicit_precomputed_quantities(Y, atmos) + FT = eltype(Y) + TST = thermo_state_type(atmos.moisture_model, FT) + n = n_mass_flux_subdomains(atmos.turbconv_model) + gs_quantities = (; + ᶜspecific = specific_gs.(Y.c), + ᶠuₕ³ = similar(Y.f, CT3{FT}), + ᶜu = similar(Y.c, C123{FT}), + ᶠu³ = similar(Y.f, CT3{FT}), + ᶜK = similar(Y.c, FT), + ᶜts = similar(Y.c, TST), + ᶜp = similar(Y.c, FT), + ᶜh_tot = similar(Y.c, FT), + ) + advective_sgs_quantities = + atmos.turbconv_model isa PrognosticEDMFX ? + (; + ᶜtke⁰ = similar(Y.c, FT), + ᶜρa⁰ = similar(Y.c, FT), + ᶜmse⁰ = similar(Y.c, FT), + ᶜq_tot⁰ = similar(Y.c, FT), + ᶠu₃⁰ = similar(Y.f, C3{FT}), + ᶜu⁰ = similar(Y.c, C123{FT}), + ᶠu³⁰ = similar(Y.f, CT3{FT}), + ᶜK⁰ = similar(Y.c, FT), + ᶜts⁰ = similar(Y.c, TST), + ᶜρ⁰ = similar(Y.c, FT), + ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}), + ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}), + ᶜKʲs = similar(Y.c, NTuple{n, FT}), + ᶠKᵥʲs = similar(Y.f, NTuple{n, FT}), + ᶜtsʲs = similar(Y.c, NTuple{n, TST}), + ᶜρʲs = similar(Y.c, NTuple{n, FT}), + ) : (;) + return (; gs_quantities..., advective_sgs_quantities...) +end + """ precomputed_quantities(Y, atmos) @@ -42,16 +79,8 @@ function precomputed_quantities(Y, atmos) TST = thermo_state_type(atmos.moisture_model, FT) SCT = SurfaceConditions.surface_conditions_type(atmos, FT) cspace = axes(Y.c) - fspace = axes(Y.f) n = n_mass_flux_subdomains(atmos.turbconv_model) gs_quantities = (; - ᶜspecific = specific_gs.(Y.c), - ᶜu = similar(Y.c, C123{FT}), - ᶠu³ = similar(Y.f, CT3{FT}), - ᶜK = similar(Y.c, FT), - ᶜts = similar(Y.c, TST), - ᶜp = similar(Y.c, FT), - ᶜh_tot = similar(Y.c, FT), ᶜmixing_length = similar(Y.c, FT), ᶜlinear_buoygrad = similar(Y.c, FT), ᶜstrain_rate_norm = similar(Y.c, FT), @@ -76,42 +105,26 @@ function precomputed_quantities(Y, atmos) advective_sgs_quantities = atmos.turbconv_model isa PrognosticEDMFX ? (; - ᶜtke⁰ = similar(Y.c, FT), - ᶜρa⁰ = similar(Y.c, FT), - ᶠu₃⁰ = similar(Y.f, C3{FT}), - ᶜu⁰ = similar(Y.c, C123{FT}), - ᶠu³⁰ = similar(Y.f, CT3{FT}), - ᶜK⁰ = similar(Y.c, FT), - ᶜmse⁰ = similar(Y.c, FT), - ᶜq_tot⁰ = similar(Y.c, FT), - ᶜts⁰ = similar(Y.c, TST), - ᶜρ⁰ = similar(Y.c, FT), ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}), ᶜK_u = similar(Y.c, FT), ᶜK_h = similar(Y.c, FT), ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}), - ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}), - ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}), - ᶜKʲs = similar(Y.c, NTuple{n, FT}), - ᶠKᵥʲs = similar(Y.f, NTuple{n, FT}), - ᶜtsʲs = similar(Y.c, NTuple{n, TST}), bdmr_l = similar(Y.c, BidiagonalMatrixRow{FT}), bdmr_r = similar(Y.c, BidiagonalMatrixRow{FT}), bdmr = similar(Y.c, BidiagonalMatrixRow{FT}), - ᶜρʲs = similar(Y.c, NTuple{n, FT}), ᶜentrʲs = similar(Y.c, NTuple{n, FT}), ᶜdetrʲs = similar(Y.c, NTuple{n, FT}), ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}), ᶠnh_pressure₃ʲs = similar(Y.f, NTuple{n, C3{FT}}), - ᶜgradᵥ_θ_virt⁰ = Fields.Field(C3{FT}, cspace), - ᶜgradᵥ_q_tot⁰ = Fields.Field(C3{FT}, cspace), - ᶜgradᵥ_θ_liq_ice⁰ = Fields.Field(C3{FT}, cspace), + ᶜgradᵥ_θ_virt⁰ = similar(Y.c, C3{FT}), + ᶜgradᵥ_q_tot⁰ = similar(Y.c, C3{FT}), + ᶜgradᵥ_θ_liq_ice⁰ = similar(Y.c, C3{FT}), precipitation_sgs_quantities..., ) : (;) sgs_quantities = (; - ᶜgradᵥ_θ_virt = Fields.Field(C3{FT}, cspace), - ᶜgradᵥ_q_tot = Fields.Field(C3{FT}, cspace), - ᶜgradᵥ_θ_liq_ice = Fields.Field(C3{FT}, cspace), + ᶜgradᵥ_θ_virt = similar(Y.c, C3{FT}), + ᶜgradᵥ_q_tot = similar(Y.c, C3{FT}), + ᶜgradᵥ_θ_liq_ice = similar(Y.c, C3{FT}), ) diagnostic_sgs_quantities = @@ -157,6 +170,7 @@ function precomputed_quantities(Y, atmos) ᶜqₛ = similar(Y.c, FT), ) : (;) return (; + implicit_precomputed_quantities(Y, atmos)..., gs_quantities..., sgs_quantities..., advective_sgs_quantities..., @@ -328,14 +342,14 @@ end function eddy_diffusivity_coefficient( z::FT, z₀, - f_b::FT, - h::FT, + f_b, + h, uₐ, - C_E::FT, - Ri::FT, - Ri_a::FT, - Ri_c::FT, - κ::FT, + C_E, + Ri, + Ri_a, + Ri_c, + κ, ) where {FT} # Equations (17), (18) if z <= f_b * h @@ -414,12 +428,12 @@ end function compute_surface_layer_diffusivity( z::FT, - z₀::FT, - κ::FT, - C_E::FT, - Ri::FT, - Ri_a::FT, - Ri_c::FT, + z₀, + κ, + C_E, + Ri, + Ri_a, + Ri_c, norm_uₐ, ) where {FT} # Equations (19), (20) @@ -454,15 +468,19 @@ function instead of recomputing the value yourself. Otherwise, it will be difficult to ensure that the duplicated computations are consistent. """ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) + set_implicit_precomputed_quantities!(Y, p, t) + set_explicit_precomputed_quantities!(Y, p, t) +end +function set_implicit_precomputed_quantities!(Y, p, t) (; moisture_model, turbconv_model, vert_diff, precip_model, cloud_model) = p.atmos (; call_cloud_diagnostics_per_stage) = p.atmos thermo_params = CAP.thermodynamics_params(p.params) + turbconv_params = CAP.turbconv_params(p.params) n = n_mass_flux_subdomains(turbconv_model) thermo_args = (thermo_params, moisture_model) (; ᶜΦ) = p.core - (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed - ᶠuₕ³ = p.scratch.ᶠtemp_CT3 + (; ᶜspecific, ᶠuₕ³, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp, ᶜh_tot) = p.precomputed @. ᶜspecific = specific_gs(Y.c) set_ᶠuₕ³!(ᶠuₕ³, Y) @@ -490,6 +508,63 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) end @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ) @. ᶜp = TD.air_pressure(thermo_params, ᶜts) + @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot) + + if turbconv_model isa PrognosticEDMFX + @assert !(moisture_model isa DryModel) + (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed + (; ᶜtke⁰, ᶜρa⁰, ᶜmse⁰, ᶜq_tot⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰) = + p.precomputed + + for j in 1:n + ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃ + ᶜmseʲ = Y.c.sgsʲs.:($j).mse + ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot + + ᶜuʲ = ᶜuʲs.:($j) + ᶠu³ʲ = ᶠu³ʲs.:($j) + ᶜKʲ = ᶜKʲs.:($j) + ᶠKᵥʲ = ᶠKᵥʲs.:($j) + ᶜtsʲ = ᶜtsʲs.:($j) + ᶜρʲ = ᶜρʲs.:($j) + + set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³) + @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2 + @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ) + @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ) + end + + @. ᶜρa⁰ = ρa⁰(Y.c) + @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model) + @. ᶜmse⁰ = divide_by_ρa( + Y.c.ρ * (ᶜh_tot - ᶜK) - ρamse⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρ * (ᶜh_tot - ᶜK), + Y.c.ρ, + turbconv_model, + ) + @. ᶜq_tot⁰ = divide_by_ρa( + Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_tot, + Y.c.ρ, + turbconv_model, + ) + set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model) + set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³) + # @. ᶜK⁰ += ᶜtke⁰ + @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰) + @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰) + end + + return nothing +end +function set_explicit_precomputed_quantities!(Y, p, t) + (; moisture_model, turbconv_model, vert_diff, precip_model, cloud_model) = + p.atmos + (; call_cloud_diagnostics_per_stage) = p.atmos + thermo_params = CAP.thermodynamics_params(p.params) + (; ᶠuₕ³, ᶜts, ᶜp) = p.precomputed if turbconv_model isa AbstractEDMF @. p.precomputed.ᶜgradᵥ_θ_virt = @@ -500,9 +575,6 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))) end - (; ᶜh_tot) = p.precomputed - @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot) - if !isnothing(p.sfc_setup) SurfaceConditions.update_surface_conditions!(Y, p, t) end @@ -519,12 +591,11 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) if turbconv_model isa PrognosticEDMFX set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t) - set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³, t) set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t) set_prognostic_edmf_precomputed_quantities_precipitation!( Y, p, - p.atmos.precip_model, + precip_model, ) end @@ -537,7 +608,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) Y, p, t, - p.atmos.precip_model, + precip_model, ) end @@ -546,7 +617,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) interior_uₕ = Fields.level(Y.c.uₕ, 1) ᶜΔz_surface = Fields.Δz_field(interior_uₕ) @. ᶜK_h = eddy_diffusivity_coefficient( - p.atmos.vert_diff.C_E, + vert_diff.C_E, norm(interior_uₕ), ᶜΔz_surface / 2, ᶜp, diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index eb620888bb..4733d53264 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -5,50 +5,6 @@ import NVTX import Thermodynamics as TD import ClimaCore: Spaces, Fields -""" - set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) - -Updates the edmf environment precomputed quantities stored in `p` for edmfx. -""" -NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!( - Y, - p, - ᶠuₕ³, - t, -) - @assert !(p.atmos.moisture_model isa DryModel) - - thermo_params = CAP.thermodynamics_params(p.params) - (; turbconv_model) = p.atmos - (; ᶜΦ,) = p.core - (; ᶜp, ᶜh_tot, ᶜK) = p.precomputed - (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜmse⁰, ᶜq_tot⁰) = - p.precomputed - - @. ᶜρa⁰ = ρa⁰(Y.c) - @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model) - @. ᶜmse⁰ = divide_by_ρa( - Y.c.ρ * (ᶜh_tot - ᶜK) - ρamse⁺(Y.c.sgsʲs), - ᶜρa⁰, - Y.c.ρ * (ᶜh_tot - ᶜK), - Y.c.ρ, - turbconv_model, - ) - @. ᶜq_tot⁰ = divide_by_ρa( - Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs), - ᶜρa⁰, - Y.c.ρq_tot, - Y.c.ρ, - turbconv_model, - ) - set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model) - set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³) - # @. ᶜK⁰ += ᶜtke⁰ - @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰) - @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰) - return nothing -end - """ set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t) @@ -62,8 +18,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! t, ) (; moisture_model, turbconv_model) = p.atmos - #EDMFX BCs only support total energy as state variable - @assert !(moisture_model isa DryModel) FT = Spaces.undertype(axes(Y.c)) n = n_mass_flux_subdomains(turbconv_model) @@ -78,21 +32,10 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions for j in 1:n - ᶜuʲ = ᶜuʲs.:($j) - ᶠu³ʲ = ᶠu³ʲs.:($j) - ᶜKʲ = ᶜKʲs.:($j) - ᶠKᵥʲ = ᶠKᵥʲs.:($j) - ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃ ᶜtsʲ = ᶜtsʲs.:($j) - ᶜρʲ = ᶜρʲs.:($j) ᶜmseʲ = Y.c.sgsʲs.:($j).mse ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot - set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³) - @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2 - @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ) - @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ) - # EDMFX boundary condition: # We need field_values everywhere because we are mixing diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl index 4482e86147..2943f52405 100644 --- a/src/cache/temporary_quantities.jl +++ b/src/cache/temporary_quantities.jl @@ -10,6 +10,7 @@ function temporary_quantities(Y, atmos) FT = Spaces.undertype(center_space) CTh = CTh_vector_type(Y.c) + num_maxtrix_rows = length(first(column_iterator(Y))) return (; ᶠtemp_scalar = Fields.Field(FT, face_space), # ᶠp, ᶠρK_E ᶜtemp_scalar = Fields.Field(FT, center_space), # ᶜ1 @@ -68,5 +69,6 @@ function temporary_quantities(Y, atmos) ᶜdiffusion_h_matrix = similar(Y.c, TridiagonalMatrixRow{FT}), ᶜdiffusion_u_matrix = similar(Y.c, TridiagonalMatrixRow{FT}), ᶠtridiagonal_matrix_c3 = similar(Y.f, TridiagonalMatrixRow{C3{FT}}), + temp_matrix = Array{FT}(undef, num_maxtrix_rows, num_maxtrix_rows), ) end diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 56bbb7076a..e9cf9b021f 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -404,7 +404,7 @@ function print_walltime_estimate(integrator) "" end estimated_finish_date = - Dates.now() + compound_period(wall_time_remaining, Dates.Second) + Dates.now() + compound_period(wall_time_remaining) @info "Progress" simulation_time = simulation_time n_steps_completed = n_steps wall_time_per_step = wall_time_ave_per_step_str wall_time_total = wall_time_total wall_time_remaining = wall_time_remaining_str wall_time_spent = @@ -452,6 +452,32 @@ function gc_func(integrator) return nothing end +function update_exact_jacobian!(integrator) + integrator.alg isa CTS.RosenbrockAlgorithm || + integrator.alg.name isa CTS.IMEXARKAlgorithmName || + return + (; tableau) = integrator.alg + tableau_coefficients = LinearAlgebra.diag( + integrator.alg isa CTS.RosenbrockAlgorithm ? tableau.Γ : tableau.a_imp, + ) + γs = unique(filter(!iszero, tableau_coefficients)) + length(γs) == 1 || + error("The exact Jacobian must be updated on every Newton iteration, \ + rather than on every timestep (or every N steps), because the \ + specified IMEX algorithm has implicit stages with distinct \ + tableau coefficients (i.e., it is not an SDIRK algorithm).") + j = + integrator.alg isa CTS.RosenbrockAlgorithm ? integrator.cache.W : + integrator.cache.newtons_method_cache.j + update_exact_jacobian!( + j, + integrator.u, + integrator.p, + integrator.dt * γs[1], + integrator.t, + ) +end + """ maybe_graceful_exit(integrator) diff --git a/src/callbacks/get_callbacks.jl b/src/callbacks/get_callbacks.jl index 45ac611530..b700d224aa 100644 --- a/src/callbacks/get_callbacks.jl +++ b/src/callbacks/get_callbacks.jl @@ -23,6 +23,7 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt) "average" => ((+), CAD.average_pre_output_hook!), ) + dict_writer = CAD.DictWriter() hdf5_writer = CAD.HDF5Writer(p.output_dir) if !isnothing(parsed_args["netcdf_interpolation_num_points"]) @@ -45,11 +46,13 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt) z_sampling_method, sync_schedule = CAD.EveryStepSchedule(), ) - writers = (hdf5_writer, netcdf_writer) - # The default writer is HDF5 + writers = (; dict_writer, hdf5_writer, netcdf_writer) + + # The default writer is NetCDF ALLOWED_WRITERS = Dict( "nothing" => netcdf_writer, + "dict" => dict_writer, "h5" => hdf5_writer, "hdf5" => hdf5_writer, "nc" => netcdf_writer, @@ -157,7 +160,30 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt) diagnostics..., ] end - diagnostics = collect(diagnostics) + if parsed_args["output_default_diagnostics"] && ( + parsed_args["use_exact_jacobian"] || + parsed_args["debug_approximate_jacobian"] + ) + period_seconds = + min( + max( + 1, + parsed_args["n_steps_update_exact_jacobian"], + fld(p.t_end / 99, p.dt), # Save no more than 100 matrices. + ), + fld(p.t_end, p.dt), # Save at least 2 matrices. + ) * p.dt + schedule = CAD.EveryDtSchedule(period_seconds) + exact_jacobian_diagnostic = CAD.ScheduledDiagnostic(; + variable = CAD.get_diagnostic_variable("ejac"), + output_schedule_func = schedule, + compute_schedule_func = deepcopy(schedule), + reduction_time_func = nothing, + pre_output_hook! = nothing, + output_writer = dict_writer, + ) + diagnostics = [exact_jacobian_diagnostic, diagnostics...] + end for writer in writers writer_str = nameof(typeof(writer)) @@ -271,6 +297,19 @@ function get_callbacks(config, sim_info, atmos, params, Y, p, t_start) (callbacks..., call_every_dt(cloud_fraction_model_callback!, dt_cf)) end + if ( + parsed_args["use_exact_jacobian"] || + parsed_args["debug_approximate_jacobian"] + ) && parsed_args["n_steps_update_exact_jacobian"] != 0 + callbacks = ( + callbacks..., + call_every_n_steps( + update_exact_jacobian!, + parsed_args["n_steps_update_exact_jacobian"], + ), + ) + end + if atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode dt_rad = FT(time_to_seconds(parsed_args["dt_rad"])) callbacks = diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 2c8a5c1928..4c8e8f008e 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -56,6 +56,7 @@ import ClimaDiagnostics.Schedules: EveryStepSchedule, EveryDtSchedule, EveryCalendarDtSchedule import ClimaDiagnostics.Writers: + DictWriter, HDF5Writer, NetCDFWriter, write_field!, diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index ab76d8a4a8..c00fba3d45 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -1163,3 +1163,19 @@ add_diagnostic_variable!( comments = "Mass of water vapor per mass of air", compute! = compute_husv!, ) + +### +# Implicit solver +### +add_diagnostic_variable!( + short_name = "ejac", + long_name = "Exact Jacobian matrix of first column", + standard_name = "exact_jacobian", + units = "s^-1", + comments = "Exact Jacobian matrix for column (1, 1, 1) of tendency", + compute! = (out, state, cache, time) -> if isnothing(out) + copy(cache.scratch.temp_matrix) + else + out .= cache.scratch.temp_matrix + end, +) diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index ac118ef08d..6880ffd090 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -108,14 +108,8 @@ end Returns the total energy source term multiplier from precipitation formation for the 0-moment scheme """ -function e_tot_0M_precipitation_sources_helper( - thp, - ts, - Φ::FT, -) where {FT <: Real} - - λ::FT = TD.liquid_fraction(thp, ts) - +function e_tot_0M_precipitation_sources_helper(thp, ts, Φ) + λ = TD.liquid_fraction(thp, ts) return λ * Iₗ(thp, ts) + (1 - λ) * Iᵢ(thp, ts) + Φ end diff --git a/src/prognostic_equations/edmfx_boundary_condition.jl b/src/prognostic_equations/edmfx_boundary_condition.jl index 20810c9df9..5d02d96799 100644 --- a/src/prognostic_equations/edmfx_boundary_condition.jl +++ b/src/prognostic_equations/edmfx_boundary_condition.jl @@ -4,9 +4,9 @@ function sgs_scalar_first_interior_bc( ᶜz_int::FT, - ᶜρ_int::FT, - ᶜaʲ_int::FT, - ᶜscalar_int::FT, + ᶜρ_int, + ᶜaʲ_int, + ᶜscalar_int, sfc_buoyancy_flux, sfc_ρ_flux_scalar, ustar, @@ -32,9 +32,9 @@ end function get_first_interior_variance( flux, - ustar::FT, + ustar, z::FT, - oblength::FT, + oblength, local_geometry, ) where {FT} c_star = -flux / ustar @@ -61,13 +61,10 @@ function guass_quantile(p::FT) where {FT <: Real} return sqrt(2) * approximate_inverf(2p - 1) end -function percentile_bounds_mean_norm( - low_percentile::FT, - high_percentile::FT, -) where {FT <: Real} +function percentile_bounds_mean_norm(low_percentile, high_percentile) gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi) xp_high = guass_quantile(high_percentile) xp_low = guass_quantile(low_percentile) return (gauss_int(xp_high) - gauss_int(xp_low)) / - max(high_percentile - low_percentile, eps(FT)) + max(high_percentile - low_percentile, eps(typeof(high_percentile))) end diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl index 018385d189..2d15310b02 100644 --- a/src/prognostic_equations/edmfx_closures.jl +++ b/src/prognostic_equations/edmfx_closures.jl @@ -153,17 +153,17 @@ Smagorinsky length scale. """ function mixing_length( params, - ustar::FT, + ustar, ᶜz::FT, - z_sfc::FT, - ᶜdz::FT, - sfc_tke::FT, - ᶜlinear_buoygrad::FT, - ᶜtke::FT, - obukhov_length::FT, - ᶜstrain_rate_norm::FT, - ᶜPr::FT, - ᶜtke_exch::FT, + z_sfc, + ᶜdz, + sfc_tke, + ᶜlinear_buoygrad, + ᶜtke, + obukhov_length, + ᶜstrain_rate_norm, + ᶜPr, + ᶜtke_exch, ) where {FT} turbconv_params = CAP.turbconv_params(params) @@ -234,7 +234,7 @@ function mixing_length( l_smin = lamb_smooth_minimum(l, smin_ub, smin_rm) l_limited = max(l_smag, min(l_smin, l_z)) - return MixingLength{FT}(l_limited, l_W, l_TKE, l_N) + return MixingLength(l_limited, l_W, l_TKE, l_N) end """ diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index 0d37e891e0..2b3148e12e 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -19,17 +19,17 @@ function entrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜtke⁰, ::NoEntrainment, ) where {FT} return FT(0) @@ -38,17 +38,17 @@ end function entrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜtke⁰, ::PiGroupsEntrainment, ) where {FT} if ᶜaʲ <= FT(0) @@ -86,17 +86,17 @@ end function entrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜtke⁰, ::GeneralizedEntrainment, ) where {FT} turbconv_params = CAP.turbconv_params(params) @@ -117,7 +117,7 @@ function entrainment( return max(entr, 0) end -function detrainment( +function detrainment_from_thermo_state( params, z_prev_level::FT, z_sfc_halflevel, @@ -177,21 +177,21 @@ end function detrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜρaʲ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜentr::FT, - ᶜvert_div::FT, - ᶜmassflux_vert_div::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜρaʲ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜentr, + ᶜvert_div, + ᶜmassflux_vert_div, + ᶜtke⁰, ::NoDetrainment, ) where {FT} return FT(0) @@ -200,21 +200,21 @@ end function detrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜρaʲ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜentr::FT, - ᶜvert_div::FT, - ᶜmassflux_vert_div::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜρaʲ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜentr, + ᶜvert_div, + ᶜmassflux_vert_div, + ᶜtke⁰, ::PiGroupsDetrainment, ) where {FT} @@ -251,21 +251,21 @@ end function detrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜρaʲ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜentr::FT, - ᶜvert_div::FT, - ᶜmassflux_vert_div::FT, - ᶜtke⁰::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜρaʲ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜentr, + ᶜvert_div, + ᶜmassflux_vert_div, + ᶜtke⁰, ::GeneralizedDetrainment, ) where {FT} turbconv_params = CAP.turbconv_params(params) @@ -301,20 +301,20 @@ end function detrainment( params, ᶜz::FT, - z_sfc::FT, - ᶜp::FT, - ᶜρ::FT, - ᶜρaʲ::FT, - ᶜaʲ::FT, - ᶜwʲ::FT, - ᶜRHʲ::FT, - ᶜbuoyʲ::FT, - ᶜw⁰::FT, - ᶜRH⁰::FT, - ᶜbuoy⁰::FT, - ᶜentr::FT, - ᶜvert_div::FT, - ᶜmassflux_vert_div::FT, + z_sfc, + ᶜp, + ᶜρ, + ᶜρaʲ, + ᶜaʲ, + ᶜwʲ, + ᶜRHʲ, + ᶜbuoyʲ, + ᶜw⁰, + ᶜRH⁰, + ᶜbuoy⁰, + ᶜentr, + ᶜvert_div, + ᶜmassflux_vert_div, ::ConstantAreaDetrainment, ) where {FT} detr = ᶜentr - ᶜvert_div diff --git a/src/prognostic_equations/gm_sgs_closures.jl b/src/prognostic_equations/gm_sgs_closures.jl index 13730270ed..c6af7fea76 100644 --- a/src/prognostic_equations/gm_sgs_closures.jl +++ b/src/prognostic_equations/gm_sgs_closures.jl @@ -19,7 +19,7 @@ function smagorinsky_lilly_length(c_smag, N_eff, dz, Pr, ϵ_st) return N_eff > FT(0) ? c_smag * dz * - max(0, 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / 4) : + max(FT(0), 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / FT(4)) : c_smag * dz end diff --git a/src/prognostic_equations/implicit/approx_jacobian.jl b/src/prognostic_equations/implicit/approx_jacobian.jl new file mode 100644 index 0000000000..079f13cd25 --- /dev/null +++ b/src/prognostic_equations/implicit/approx_jacobian.jl @@ -0,0 +1,683 @@ +abstract type DerivativeFlag end +struct UseDerivative <: DerivativeFlag end +struct IgnoreDerivative <: DerivativeFlag end +use_derivative(::UseDerivative) = true +use_derivative(::IgnoreDerivative) = false + +""" + ApproxJacobian(Y, p; approximate_solve_iters, diffusion_flag, topography_flag) + +An approximation of the exact `ImplicitEquationJacobian`, which is updated using +analytically derived tendency derivatives and inverted using a specialized +nested linear solver. + +# Keyword Arguments + +- `diffusion_flag::DerivativeFlag`: whether the derivative of the + diffusion tendency with respect to the quantities being diffused should be + computed or approximated as 0; must be either `UseDerivative()` or + `Ignoreerivative()` instead of a `Bool` to ensure type-stability +- `topography_flag::DerivativeFlag`: whether the derivative of vertical + contravariant velocity with respect to horizontal covariant velocity should be + computed or approximated as 0; must be either `UseDerivative()` or + `IgnoreDerivative()` instead of a `Bool` to ensure type-stability +- `sgs_advection_flag::DerivativeFlag`: whether the derivative of the + subgrid-scale advection tendency with respect to the updraft quantities should + be computed or approximated as 0; must be either `UseDerivative()` or + `IgnoreDerivative()` instead of a `Bool` to ensure type-stability +- `approximate_solve_iters::Int`: number of iterations to take for the + approximate linear solve required when `diffusion_flag = UseDerivative()` +""" +@kwdef struct ApproxJacobian{F1, F2, F3} <: JacobianAlgorithm + diffusion_flag::F1 = nothing + topography_flag::F2 = nothing + sgs_advection_flag::F3 = nothing + approximate_solve_iters::Int = 1 +end + +function jacobian_cache(alg::ApproxJacobian, Y, p) + diffusion_flag = if isnothing(alg.diffusion_flag) + p.atmos.diff_mode == Implicit() ? UseDerivative() : IgnoreDerivative() + else + alg.diffusion_flag + end + topography_flag = if isnothing(alg.topography_flag) + has_topography(axes(Y.c)) ? UseDerivative() : IgnoreDerivative() + else + alg.topography_flag + end + sgs_advection_flag = if isnothing(alg.sgs_advection_flag) + p.atmos.sgs_adv_mode == Implicit() ? UseDerivative() : IgnoreDerivative() + else + alg.sgs_advection_flag + end + + FT = Spaces.undertype(axes(Y.c)) + CTh = CTh_vector_type(axes(Y.c)) + + DiagonalRow = DiagonalMatrixRow{FT} + TridiagonalRow = TridiagonalMatrixRow{FT} + BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}} + TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}} + BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}} + BidiagonalRow_C3xACTh = + BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CTh{FT})')} + DiagonalRow_C3xACT3 = + DiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} + TridiagonalRow_C3xACT3 = + TridiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} + + is_in_Y(name) = MatrixFields.has_field(Y, name) + + ρq_tot_if_available = is_in_Y(@name(c.ρq_tot)) ? (@name(c.ρq_tot),) : () + ρatke_if_available = + is_in_Y(@name(c.sgs⁰.ρatke)) ? (@name(c.sgs⁰.ρatke),) : () + sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : () + + tracer_names = ( + @name(c.ρq_tot), + @name(c.ρq_liq), + @name(c.ρq_ice), + @name(c.ρq_rai), + @name(c.ρq_sno), + ) + available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names) + + # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0, + # which means that multiplying inv(-1) by a Float32 will yield a Float64. + identity_blocks = MatrixFields.unrolled_map( + name -> (name, name) => FT(-1) * I, + (@name(c.ρ), sfc_if_available...), + ) + + active_scalar_names = (@name(c.ρ), @name(c.ρe_tot), ρq_tot_if_available...) + advection_blocks = ( + ( + use_derivative(topography_flag) ? + MatrixFields.unrolled_map( + name -> + (name, @name(c.uₕ)) => + similar(Y.c, TridiagonalRow_ACTh), + active_scalar_names, + ) : () + )..., + MatrixFields.unrolled_map( + name -> (name, @name(f.u₃)) => similar(Y.c, BidiagonalRow_ACT3), + active_scalar_names, + )..., + MatrixFields.unrolled_map( + name -> (@name(f.u₃), name) => similar(Y.f, BidiagonalRow_C3), + active_scalar_names, + )..., + (@name(f.u₃), @name(c.uₕ)) => similar(Y.f, BidiagonalRow_C3xACTh), + (@name(f.u₃), @name(f.u₃)) => similar(Y.f, TridiagonalRow_C3xACT3), + ) + + diffused_scalar_names = + (@name(c.ρe_tot), available_tracer_names..., ρatke_if_available...) + diffusion_blocks = if use_derivative(diffusion_flag) + ( + MatrixFields.unrolled_map( + name -> (name, @name(c.ρ)) => similar(Y.c, TridiagonalRow), + diffused_scalar_names, + )..., + MatrixFields.unrolled_map( + name -> (name, name) => similar(Y.c, TridiagonalRow), + diffused_scalar_names, + )..., + ( + is_in_Y(@name(c.ρq_tot)) ? + ( + (@name(c.ρe_tot), @name(c.ρq_tot)) => + similar(Y.c, TridiagonalRow), + ) : () + )..., + (@name(c.uₕ), @name(c.uₕ)) => + !isnothing(p.atmos.turbconv_model) || + diffuse_momentum(p.atmos.vert_diff) ? + similar(Y.c, TridiagonalRow) : FT(-1) * I, + ) + else + MatrixFields.unrolled_map( + name -> (name, name) => FT(-1) * I, + (diffused_scalar_names..., @name(c.uₕ)), + ) + end + + sgs_advection_blocks = if p.atmos.turbconv_model isa PrognosticEDMFX + @assert n_prognostic_mass_flux_subdomains(p.atmos.turbconv_model) == 1 + sgs_scalar_names = ( + @name(c.sgsʲs.:(1).q_tot), + @name(c.sgsʲs.:(1).mse), + @name(c.sgsʲs.:(1).ρa) + ) + if use_derivative(sgs_advection_flag) + ( + MatrixFields.unrolled_map( + name -> (name, name) => similar(Y.c, TridiagonalRow), + sgs_scalar_names, + )..., + (@name(c.sgsʲs.:(1).mse), @name(c.ρ)) => + similar(Y.c, DiagonalRow), + (@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.c, DiagonalRow), + (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.c, TridiagonalRow), + (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)) => + similar(Y.c, TridiagonalRow), + (@name(f.sgsʲs.:(1).u₃), @name(c.ρ)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => + similar(Y.f, TridiagonalRow_C3xACT3), + ) + else + ( + MatrixFields.unrolled_map( + name -> (name, name) => FT(-1) * I, + sgs_scalar_names, + )..., + (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => + !isnothing(p.atmos.rayleigh_sponge) ? + similar(Y.f, DiagonalRow_C3xACT3) : FT(-1) * I, + ) + end + else + () + end + + matrix = MatrixFields.FieldMatrix( + identity_blocks..., + sgs_advection_blocks..., + advection_blocks..., + diffusion_blocks..., + ) + + sgs_names_if_available = if p.atmos.turbconv_model isa PrognosticEDMFX + ( + @name(c.sgsʲs.:(1).q_tot), + @name(c.sgsʲs.:(1).mse), + @name(c.sgsʲs.:(1).ρa), + @name(f.sgsʲs.:(1).u₃), + ) + else + () + end + names₁_group₁ = (@name(c.ρ), sfc_if_available...) + names₁_group₂ = (available_tracer_names..., ρatke_if_available...) + names₁_group₃ = (@name(c.ρe_tot),) + names₁ = ( + names₁_group₁..., + sgs_names_if_available..., + names₁_group₂..., + names₁_group₃..., + ) + + alg₂ = MatrixFields.BlockLowerTriangularSolve(@name(c.uₕ)) + alg = + if use_derivative(diffusion_flag) || use_derivative(sgs_advection_flag) + alg₁_subalg₂ = + if p.atmos.turbconv_model isa PrognosticEDMFX && + use_derivative(sgs_advection_flag) + diff_subalg = + use_derivative(diffusion_flag) ? + (; + alg₂ = MatrixFields.BlockLowerTriangularSolve( + names₁_group₂..., + ) + ) : (;) + (; + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).q_tot); + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).mse); + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).ρa), + @name(f.sgsʲs.:(1).u₃); + diff_subalg..., + ), + ), + ) + ) + else + is_in_Y(@name(c.ρq_tot)) ? + (; + alg₂ = MatrixFields.BlockLowerTriangularSolve( + names₁_group₂..., + ) + ) : (;) + end + alg₁ = MatrixFields.BlockLowerTriangularSolve( + names₁_group₁...; + alg₁_subalg₂..., + ) + MatrixFields.ApproximateBlockArrowheadIterativeSolve( + names₁...; + alg₁, + alg₂, + P_alg₁ = MatrixFields.MainDiagonalPreconditioner(), + n_iters = approximate_solve_iters, + ) + else + MatrixFields.BlockArrowheadSolve(names₁...; alg₂) + end + + solver = MatrixFields.FieldMatrixSolver(alg, matrix, Y) + + return (; + diffusion_flag, + topography_flag, + sgs_advection_flag, + matrix, + solver, + ) +end + +always_update_exact_jacobian(::ApproxJacobian) = false + +factorize_exact_jacobian!(::ApproxJacobian, _, _, _, _, _) = nothing + +function approximate_jacobian!(::ApproxJacobian, cache, Y, p, dtγ, t) + (; diffusion_flag, topography_flag, sgs_advection_flag, matrix) = cache + (; params) = p + (; edmfx_upwinding) = p.atmos.numerics + (; ᶜΦ, ᶠgradᵥ_ᶜΦ) = p.core + (; ᶜspecific, ᶜK, ᶜts, ᶜp) = p.precomputed + (; + ∂ᶜK_∂ᶜuₕ, + ∂ᶜK_∂ᶠu₃, + ᶠp_grad_matrix, + ᶜadvection_matrix, + ᶜdiffusion_h_matrix, + ᶜdiffusion_u_matrix, + ᶠbidiagonal_matrix_ct3, + ᶠbidiagonal_matrix_ct3_2, + ᶠtridiagonal_matrix_c3, + ) = p.scratch + + FT = Spaces.undertype(axes(Y.c)) + CTh = CTh_vector_type(axes(Y.c)) + one_C3xACT3 = C3(FT(1)) * CT3(FT(1))' + + cv_d = FT(CAP.cv_d(params)) + Δcv_v = FT(CAP.cv_v(params)) - cv_d + T_tri = FT(CAP.T_triple(params)) + e_int_v0 = T_tri * Δcv_v - FT(CAP.e_int_v0(params)) + thermo_params = CAP.thermodynamics_params(params) + + ᶜρ = Y.c.ρ + ᶜuₕ = Y.c.uₕ + ᶠu₃ = Y.f.u₃ + ᶜJ = Fields.local_geometry_field(Y.c).J + ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ + ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ + + ᶜkappa_m = p.scratch.ᶜtemp_scalar + @. ᶜkappa_m = + TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts) + + if use_derivative(topography_flag) + @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow( + adjoint(CTh(ᶜuₕ)) + adjoint(ᶜinterp(ᶠu₃)) * g³ʰ(ᶜgⁱʲ), + ) + else + @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ))) + end + @. ∂ᶜK_∂ᶠu₃ = + ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) + + DiagonalMatrixRow(adjoint(CT3(ᶜuₕ))) ⋅ ᶜinterp_matrix() + + @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() + + @. ᶜadvection_matrix = + -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) + + if use_derivative(topography_flag) + ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)] + @. ∂ᶜρ_err_∂ᶜuₕ = + dtγ * ᶜadvection_matrix ⋅ ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ + DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ)) + end + ∂ᶜρ_err_∂ᶠu₃ = matrix[@name(c.ρ), @name(f.u₃)] + @. ∂ᶜρ_err_∂ᶠu₃ = dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ)) + + tracer_info = ( + (@name(c.ρe_tot), @name(ᶜh_tot)), + (@name(c.ρq_tot), @name(ᶜspecific.q_tot)), + ) + MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name) + MatrixFields.has_field(Y, ρχ_name) || return + ᶜχ = MatrixFields.get_field(p.precomputed, χ_name) + if use_derivative(topography_flag) + ∂ᶜρχ_err_∂ᶜuₕ = matrix[ρχ_name, @name(c.uₕ)] + end + ∂ᶜρχ_err_∂ᶠu₃ = matrix[ρχ_name, @name(f.u₃)] + use_derivative(topography_flag) && @. ∂ᶜρχ_err_∂ᶜuₕ = + dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅ + ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ)) + @. ∂ᶜρχ_err_∂ᶠu₃ = + dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ) * g³³(ᶠgⁱʲ)) + end + + ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)] + ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)] + @. ∂ᶠu₃_err_∂ᶜρ = + dtγ * ( + ᶠp_grad_matrix ⋅ + DiagonalMatrixRow(ᶜkappa_m * (T_tri * cv_d - ᶜK - ᶜΦ)) + + DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅ + ᶠinterp_matrix() + ) + @. ∂ᶠu₃_err_∂ᶜρe_tot = dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m) + if MatrixFields.has_field(Y, @name(c.ρq_tot)) + ∂ᶠu₃_err_∂ᶜρq_tot = matrix[@name(f.u₃), @name(c.ρq_tot)] + @. ∂ᶠu₃_err_∂ᶜρq_tot = + dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m * e_int_v0) + end + + ∂ᶠu₃_err_∂ᶜuₕ = matrix[@name(f.u₃), @name(c.uₕ)] + ∂ᶠu₃_err_∂ᶠu₃ = matrix[@name(f.u₃), @name(f.u₃)] + I_u₃ = DiagonalMatrixRow(one_C3xACT3) + @. ∂ᶠu₃_err_∂ᶜuₕ = + dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ ∂ᶜK_∂ᶜuₕ + if p.atmos.rayleigh_sponge isa RayleighSponge + @. ∂ᶠu₃_err_∂ᶠu₃ = + dtγ * ( + ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ + ∂ᶜK_∂ᶠu₃ + DiagonalMatrixRow(-p.ᶠβ_rayleigh_w * (one_C3xACT3,)) + ) - (I_u₃,) + else + @. ∂ᶠu₃_err_∂ᶠu₃ = + dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ + ∂ᶜK_∂ᶠu₃ - (I_u₃,) + end + + if use_derivative(diffusion_flag) + (; ᶜK_h, ᶜK_u) = p + @. ᶜdiffusion_h_matrix = + ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅ + ᶠgradᵥ_matrix() + if ( + MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) || + !isnothing(p.atmos.turbconv_model) || + diffuse_momentum(p.atmos.vert_diff) + ) + @. ᶜdiffusion_u_matrix = + ᶜadvdivᵥ_matrix() ⋅ + DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)) ⋅ ᶠgradᵥ_matrix() + end + + ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name(c.ρe_tot), @name(c.ρ)] + ∂ᶜρe_tot_err_∂ᶜρe_tot = matrix[@name(c.ρe_tot), @name(c.ρe_tot)] + @. ∂ᶜρe_tot_err_∂ᶜρ = + dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow( + ( + -(1 + ᶜkappa_m) * ᶜspecific.e_tot - + ᶜkappa_m * e_int_v0 * ᶜspecific.q_tot + ) / ᶜρ, + ) + @. ∂ᶜρe_tot_err_∂ᶜρe_tot = + dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow((1 + ᶜkappa_m) / ᶜρ) - + (I,) + if MatrixFields.has_field(Y, @name(c.ρq_tot)) + ∂ᶜρe_tot_err_∂ᶜρq_tot = matrix[@name(c.ρe_tot), @name(c.ρq_tot)] + @. ∂ᶜρe_tot_err_∂ᶜρq_tot = + dtγ * ᶜdiffusion_h_matrix ⋅ + DiagonalMatrixRow(ᶜkappa_m * e_int_v0 / ᶜρ) + end + + tracer_info = ( + (@name(c.ρq_tot), @name(q_tot)), + (@name(c.ρq_liq), @name(q_liq)), + (@name(c.ρq_ice), @name(q_ice)), + (@name(c.ρq_rai), @name(q_rai)), + (@name(c.ρq_sno), @name(q_sno)), + ) + MatrixFields.unrolled_foreach(tracer_info) do (ρq_name, q_name) + MatrixFields.has_field(Y, ρq_name) || return + ᶜq = MatrixFields.get_field(ᶜspecific, q_name) + ∂ᶜρq_err_∂ᶜρ = matrix[ρq_name, @name(c.ρ)] + ∂ᶜρq_err_∂ᶜρq = matrix[ρq_name, ρq_name] + @. ∂ᶜρq_err_∂ᶜρ = + dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ) + @. ∂ᶜρq_err_∂ᶜρq = + dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ) - (I,) + end + + if MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) + turbconv_params = CAP.turbconv_params(params) + c_d = CAP.tke_diss_coeff(turbconv_params) + (; ᶜtke⁰, ᶜmixing_length, dt) = p + ᶜρa⁰ = p.atmos.turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : ᶜρ + ᶜρatke⁰ = Y.c.sgs⁰.ρatke + + @inline dissipation_rate(tke⁰, mixing_length) = + tke⁰ >= 0 ? c_d * sqrt(tke⁰) / max(mixing_length, 1) : 1 / dt + @inline ∂dissipation_rate_∂tke⁰(tke⁰, mixing_length) = + tke⁰ > 0 ? c_d / (2 * max(mixing_length, 1) * sqrt(tke⁰)) : + typeof(tke⁰)(0) + + ᶜdissipation_matrix_diagonal = p.scratch.ᶜtemp_scalar + @. ᶜdissipation_matrix_diagonal = + ᶜρatke⁰ * ∂dissipation_rate_∂tke⁰(ᶜtke⁰, ᶜmixing_length) + + ∂ᶜρatke⁰_err_∂ᶜρ = matrix[@name(c.sgs⁰.ρatke), @name(c.ρ)] + ∂ᶜρatke⁰_err_∂ᶜρatke⁰ = + matrix[@name(c.sgs⁰.ρatke), @name(c.sgs⁰.ρatke)] + @. ∂ᶜρatke⁰_err_∂ᶜρ = + dtγ * ( + ᶜdiffusion_u_matrix - + DiagonalMatrixRow(ᶜdissipation_matrix_diagonal) + ) ⋅ DiagonalMatrixRow(-(ᶜtke⁰) / ᶜρa⁰) + @. ∂ᶜρatke⁰_err_∂ᶜρatke⁰ = + dtγ * ( + ( + ᶜdiffusion_u_matrix - + DiagonalMatrixRow(ᶜdissipation_matrix_diagonal) + ) ⋅ DiagonalMatrixRow(1 / ᶜρa⁰) - + DiagonalMatrixRow(dissipation_rate(ᶜtke⁰, ᶜmixing_length)) + ) - (I,) + end + + if ( + !isnothing(p.atmos.turbconv_model) || + diffuse_momentum(p.atmos.vert_diff) + ) + ∂ᶜuₕ_err_∂ᶜuₕ = matrix[@name(c.uₕ), @name(c.uₕ)] + @. ∂ᶜuₕ_err_∂ᶜuₕ = + dtγ * DiagonalMatrixRow(1 / ᶜρ) ⋅ ᶜdiffusion_u_matrix - (I,) + end + + ᶠlg = Fields.local_geometry_field(Y.f) + precip_info = + ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ))) + MatrixFields.unrolled_foreach(precip_info) do (ρqₚ_name, wₚ_name) + MatrixFields.has_field(Y, ρqₚ_name) || return + ∂ᶜρqₚ_err_∂ᶜρqₚ = matrix[ρqₚ_name, ρqₚ_name] + ᶜwₚ = MatrixFields.get_field(p.precomputed, wₚ_name) + ᶠtmp = p.scratch.ᶠtemp_CT3 + @. ᶠtmp = CT3(unit_basis_vector_data(CT3, ᶠlg)) * ᶠwinterp(ᶜJ, ᶜρ) + @. ∂ᶜρqₚ_err_∂ᶜρqₚ += + dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅ + ᶠright_bias_matrix() ⋅ DiagonalMatrixRow(-(ᶜwₚ) / ᶜρ) + end + end + + if p.atmos.turbconv_model isa PrognosticEDMFX + if use_derivative(sgs_advection_flag) + (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p + (; bdmr_l, bdmr_r, bdmr) = p + is_third_order = edmfx_upwinding == Val(:third_order) + ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1 + ᶠset_upwind_bcs = Operators.SetBoundaryOperator(; + top = Operators.SetValue(zero(CT3{FT})), + bottom = Operators.SetValue(zero(CT3{FT})), + ) # Need to wrap ᶠupwind in this for well-defined boundaries. + UpwindMatrixRowType = + is_third_order ? QuaddiagonalMatrixRow : BidiagonalMatrixRow + ᶠupwind_matrix = is_third_order ? ᶠupwind3_matrix : ᶠupwind1_matrix + ᶠset_upwind_matrix_bcs = Operators.SetBoundaryOperator(; + top = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), + bottom = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), + ) # Need to wrap ᶠupwind_matrix in this for well-defined boundaries. + ᶜkappa_mʲ = p.scratch.ᶜtemp_scalar + @. ᶜkappa_mʲ = + TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)) / + TD.cv_m(thermo_params, ᶜtsʲs.:(1)) + ∂ᶜq_totʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ = + dtγ * ( + DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) - + ᶜadvdivᵥ_matrix() ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) + ) - (I,) + + ∂ᶜmseʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶜmseʲ_err_∂ᶜq_totʲ = + dtγ * ( + -DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * + ᶜgradᵥ_ᶠΦ * + Y.c.ρ * + ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * e_int_v0, + ) + ) + ∂ᶜmseʲ_err_∂ᶜρ = matrix[@name(c.sgsʲs.:(1).mse), @name(c.ρ)] + @. ∂ᶜmseʲ_err_∂ᶜρ = + dtγ * ( + -DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * ᶜgradᵥ_ᶠΦ / ᶜρʲs.:(1), + ) + ) + ∂ᶜmseʲ_err_∂ᶜmseʲ = + matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).mse)] + @. ∂ᶜmseʲ_err_∂ᶜmseʲ = + dtγ * ( + DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) - + ᶜadvdivᵥ_matrix() ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) - + DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * + ᶜgradᵥ_ᶠΦ * + Y.c.ρ * + ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp), + ) + ) - (I,) + + ∂ᶜρaʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)] + @. ᶠbidiagonal_matrix_ct3 = + DiagonalMatrixRow( + ᶠset_upwind_bcs( + ᶠupwind( + ᶠu³ʲs.:(1), + draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)), + ), + ), + ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) * + e_int_v0, + ) + @. ᶠbidiagonal_matrix_ct3_2 = + DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ + DiagonalMatrixRow( + Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * + e_int_v0, + ) + + @. ∂ᶜρaʲ_err_∂ᶜq_totʲ = + dtγ * ᶜadvdivᵥ_matrix() ⋅ + (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2) + ∂ᶜρaʲ_err_∂ᶜmseʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)] + @. ᶠbidiagonal_matrix_ct3 = + DiagonalMatrixRow( + ᶠset_upwind_bcs( + ᶠupwind( + ᶠu³ʲs.:(1), + draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)), + ), + ), + ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp), + ) + @. ᶠbidiagonal_matrix_ct3_2 = + DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ + DiagonalMatrixRow( + Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp), + ) + @. ∂ᶜρaʲ_err_∂ᶜmseʲ = + dtγ * ᶜadvdivᵥ_matrix() ⋅ + (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2) + ∂ᶜρaʲ_err_∂ᶜρaʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)] + @. ᶜadvection_matrix = + -(ᶜadvdivᵥ_matrix()) ⋅ + DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) + @. ∂ᶜρaʲ_err_∂ᶜρaʲ = + dtγ * ᶜadvection_matrix ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ + DiagonalMatrixRow(1 / ᶜρʲs.:(1)) - (I,) + + ∂ᶠu₃ʲ_err_∂ᶜρ = matrix[@name(f.sgsʲs.:(1).u₃), @name(c.ρ)] + @. ∂ᶠu₃ʲ_err_∂ᶜρ = + dtγ * DiagonalMatrixRow(ᶠgradᵥ_ᶜΦ / ᶠinterp(ᶜρʲs.:(1))) ⋅ + ᶠinterp_matrix() + ∂ᶠu₃ʲ_err_∂ᶜq_totʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ = + dtγ * DiagonalMatrixRow( + ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2, + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) * + e_int_v0, + ) + ∂ᶠu₃ʲ_err_∂ᶜmseʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)] + @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ = + dtγ * DiagonalMatrixRow( + ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2, + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp), + ) + ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] + ᶜu₃ʲ = p.scratch.ᶜtemp_C3 + @. ᶜu₃ʲ = ᶜinterp(Y.f.sgsʲs.:(1).u₃) + + @. bdmr_l = convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix()) + @. bdmr_r = convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix()) + @. bdmr = ifelse(ᶜu₃ʲ.components.data.:1 > 0, bdmr_l, bdmr_r) + + @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋅ bdmr + if p.atmos.rayleigh_sponge isa RayleighSponge + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + dtγ * ( + ᶠtridiagonal_matrix_c3 ⋅ + DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - + DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) + ) - (I_u₃,) + else + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + dtγ * ᶠtridiagonal_matrix_c3 ⋅ + DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - (I_u₃,) + end + elseif p.atmos.rayleigh_sponge isa RayleighSponge + ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + dtγ * -DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) - + (I_u₃,) + end + end +end + +invert_jacobian!(::ApproxJacobian, cache, x, b) = + MatrixFields.field_matrix_solve!(cache.solver, x, cache.matrix, b) diff --git a/src/prognostic_equations/implicit/debug_jacobian.jl b/src/prognostic_equations/implicit/debug_jacobian.jl new file mode 100644 index 0000000000..4be4a13018 --- /dev/null +++ b/src/prognostic_equations/implicit/debug_jacobian.jl @@ -0,0 +1,108 @@ +struct DebugJacobian{E <: ExactJacobian, A <: ApproxJacobian} <: + JacobianAlgorithm + exact_jacobian_mode::E + approx_jacobian_mode::A + use_exact_jacobian::Bool +end +function DebugJacobian( + approx_jacobian_mode; + use_exact_jacobian = true, + exact_jacobian_mode_kwargs..., +) + exact_jacobian_mode = ExactJacobian(; + preserve_unfactorized_jacobian = true, + exact_jacobian_mode_kwargs..., + ) + return DebugJacobian( + exact_jacobian_mode, + approx_jacobian_mode, + use_exact_jacobian, + ) +end + +jacobian_cache(alg::DebugJacobian, Y, p) = (; + exact_jacobian_cache = jacobian_cache(alg.exact_jacobian_mode, Y, p), + approx_jacobian_cache = jacobian_cache(alg.approx_jacobian_mode, Y, p), + similar_to_x = similar(Y), + x_error = similar(Y), + x_rel_error = similar(Y), + compute_error_ref = Ref{Bool}(), + t_ref = Ref{typeof(p.t_end)}(), +) + +always_update_exact_jacobian(alg::DebugJacobian) = + always_update_exact_jacobian(alg.exact_jacobian_mode) + +function factorize_exact_jacobian!(alg::DebugJacobian, cache, Y, p, dtγ, t) + factorize_exact_jacobian!( + alg.exact_jacobian_mode, + cache.exact_jacobian_cache, + Y, + p, + dtγ, + t, + ) + cache.compute_error_ref[] = true + cache.t_ref[] = t +end + +approximate_jacobian!(alg::DebugJacobian, cache, Y, p, dtγ, t) = + (!alg.use_exact_jacobian || cache.compute_error_ref[]) && + approximate_jacobian!( + alg.approx_jacobian_mode, + cache.approx_jacobian_cache, + Y, + p, + dtγ, + t, + ) + +import Dates: @sprintf # TODO: Add Printf as a dependency of ClimaAtmos. + +function invert_jacobian!(alg::DebugJacobian, cache, x, b) + (; exact_jacobian_mode, approx_jacobian_mode, use_exact_jacobian) = alg + (; exact_jacobian_cache, approx_jacobian_cache) = cache + (; similar_to_x, x_error, x_rel_error, compute_error_ref, t_ref) = cache + if compute_error_ref[] + x_exact = use_exact_jacobian ? x : similar_to_x + x_approx = use_exact_jacobian ? similar_to_x : x + x_exact .= x_approx .= NaN # TODO: Remove this sanity check + invert_jacobian!(exact_jacobian_mode, exact_jacobian_cache, x_exact, b) + invert_jacobian!( + approx_jacobian_mode, + approx_jacobian_cache, + x_approx, + b, + ) + compute_error_ref[] = false + + @. x_error = abs(x_approx - x_exact) + @. x_rel_error = ifelse( + x_approx == x_exact == 0, + 0, + abs((x_approx - x_exact) / x_exact), + ) + + rms(f, x) = sqrt(mean(abs2 ∘ f, x)) + rms(x) = rms(identity, x) + + field_vector_rms(x) = rms(rms ∘ Fields.backing_array, Fields._values(x)) + field_vector_max(x) = + maximum(maximum ∘ Fields.backing_array, Fields._values(x)) + + rel_rms_error = field_vector_rms(x_error) / field_vector_rms(x_exact) + rel_max_error = field_vector_max(x_error) / field_vector_max(x_exact) + rms_rel_error = field_vector_rms(x_rel_error) + max_rel_error = field_vector_max(x_rel_error) + @info "Error of approximate implicit solve at t = $(t_ref[]):\n \ + relative RMS error = $(@sprintf("%.3e", rel_rms_error))\n \ + relative max error = $(@sprintf("%.3e", rel_max_error))\n \ + RMS relative error = $(@sprintf("%.3e", rms_rel_error))\n \ + max relative error = $(@sprintf("%.3e", max_rel_error))" + @assert rel_rms_error < 1 # TODO: Make this test more rigorous + elseif use_exact_jacobian + invert_jacobian!(exact_jacobian_mode, exact_jacobian_cache, x, b) + else + invert_jacobian!(approx_jacobian_mode, approx_jacobian_cache, x, b) + end +end diff --git a/src/prognostic_equations/implicit/dual_fixes.jl b/src/prognostic_equations/implicit/dual_fixes.jl new file mode 100644 index 0000000000..2dffd03289 --- /dev/null +++ b/src/prognostic_equations/implicit/dual_fixes.jl @@ -0,0 +1,182 @@ +# TODO: Move these to Thermodynamics.jl +TD.air_density(param_set, T, p, q) = + TD.air_density(param_set, TD.promote_phase_partition(T, p, q)...) +TD.air_pressure(param_set, T, ρ, q) = + TD.air_pressure(param_set, TD.promote_phase_partition(T, ρ, q)...) +TD.internal_energy(param_set, T, q) = + TD.internal_energy(param_set, TD.promote_phase_partition(T, T, q)[2:end]...) +TD.internal_energy_sat(param_set, T, ρ, q_tot, type) = + TD.internal_energy_sat(param_set, promote(T, ρ, q_tot)..., type) +TD.PhasePartition_equil(param_set, T, ρ, q_tot, type) = + TD.PhasePartition_equil(param_set, promote(T, ρ, q_tot)..., type) +TD.PhasePartition_equil_given_p(param_set, T, p, q_tot, type) = + TD.PhasePartition_equil_given_p(param_set, promote(T, p, q_tot)..., type) +TD.specific_enthalpy_sat(param_set, T, ρ, q_tot, type) = + TD.specific_enthalpy_sat(param_set, promote(T, ρ, q_tot)..., type) +TD.saturation_vapor_pressure(param_set, T, LH_0, Δcp) = + TD.saturation_vapor_pressure(param_set, promote(T, LH_0, Δcp)...) + +# TODO: Move these to SurfaceFluxes.jl +import SurfaceFluxes as SF +SF.FluxesAndFrictionVelocity(int_values, surf_values, args...) = + SF.FluxesAndFrictionVelocity(int_values, surf_values, promote(args...)...) + +# TODO: Move this to ClimaCore.jl +ClimaComms.device(fv::Fields.FieldVector) = + ClimaComms.device(first(Fields._values(fv))) + +# Ensure that we can broadcast between FieldVectors and other BlockVectors, but +# without going through the type-unstable BlockArrays broadcasting code. +# TODO: Replace the FieldVector broadcasting code in ClimaCore.jl with this +struct NewFieldVectorStyle <: Base.Broadcast.AbstractArrayStyle{1} end +Base.Broadcast.BroadcastStyle(::Type{<:Fields.FieldVector}) = + NewFieldVectorStyle() +Base.Broadcast.BroadcastStyle( + fs::NewFieldVectorStyle, + as::Base.Broadcast.DefaultArrayStyle{N}, +) where {N} = N in (0, 1) ? fs : as # prevents an ambiguity with Base.Broadcast +Base.Broadcast.BroadcastStyle( + fs::NewFieldVectorStyle, + as::Base.Broadcast.AbstractArrayStyle{N}, +) where {N} = N in (0, 1) ? fs : as +Base.@propagate_inbounds function BlockArrays.viewblock( + fv::Fields.FieldVector, + block::BlockArrays.Block{1}, +) + array = Fields.backing_array(Fields._values(fv)[block.n...]) + # vec(array) allocates; see https://github.com/JuliaLang/julia/issues/36313 + return Base.ReshapedArray(array, (length(array),), ()) +end + +first_field_vector(x::Fields.FieldVector, args...) = x +@inline first_field_vector(x::Base.Broadcast.Broadcasted, args...) = + first_field_vector(x.args..., args...) +@inline first_field_vector(_, args...) = first_field_vector(args...) +check_broadcast_axes(fv, x::Base.Broadcast.Broadcasted) = + check_broadcast_args_axes(fv, x.args...) +@inline check_broadcast_axes(fv, x) = + isempty(axes(x)) || # scalar + axes(x) == axes(fv) || # AbstractBlockVector + x isa AbstractVector && length(x) == length(fv) || # fallback + error("Mismatched broadcast axes: $(axes(x)) != $(axes(fv))") +check_broadcast_args_axes(_) = true +@inline check_broadcast_args_axes(fv, arg, args...) = + check_broadcast_axes(fv, arg) && check_broadcast_args_axes(fv, args...) + +block_view(x, _, _) = x +Base.@propagate_inbounds block_view( + x::BlockArrays.AbstractBlockVector, + ::Val{index}, + _, +) where {index} = BlockArrays.viewblock(x, BlockArrays.Block(index)) +Base.@propagate_inbounds block_view(x::AbstractVector, val, fv) = + block_view(BlockArrays.PseudoBlockVector(x, axes(fv)), val, fv) +Base.@propagate_inbounds block_view(x::Base.Broadcast.Broadcasted, val, fv) = + Base.Broadcast.broadcasted(x.f, block_view_args(val, fv, x.args...)...) +block_view_args(_, _) = () +Base.@propagate_inbounds block_view_args(val, fv, arg, args...) = + (block_view(arg, val, fv), block_view_args(val, fv, args...)...) + +# Bypass the method for materialize! defined in BlockArrays, which results in +# type instabilities and allocations. +Base.materialize!( + dest::Fields.FieldVector, + bc::Base.Broadcast.Broadcasted{NewFieldVectorStyle}, +) = materialize_field_vector!(dest, bc) +Base.materialize!( + dest::AbstractVector, + bc::Base.Broadcast.Broadcasted{NewFieldVectorStyle}, +) = materialize_field_vector!(dest, bc) +Base.materialize!( + dest::Fields.FieldVector, + bc::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{0}}, +) = materialize_field_vector!(dest, bc) +Base.materialize!( + dest::Fields.FieldVector, + bc::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{1}}, +) = materialize_field_vector!(dest, bc) +Base.materialize!( + dest::Fields.FieldVector, + bc::Base.Broadcast.Broadcasted{<:BlockArrays.AbstractBlockStyle{1}}, +) = materialize_field_vector!(dest, bc) # prevents an ambiguity with BlockArrays +function materialize_field_vector!(dest, bc) + fv = first_field_vector(dest, bc) + check_broadcast_args_axes(fv, dest, bc) + block_index_vals = ntuple(Val, BlockArrays.blocksize(fv, 1)) + MatrixFields.unrolled_foreach(block_index_vals) do val + dest_view = block_view(dest, val, fv) + bc_view = block_view(bc, val, fv) + @inbounds Base.materialize!(dest_view, bc_view) + end + return dest +end +Base.@propagate_inbounds function Base.getindex( + fv::Fields.FieldVector, + i::Integer, +) + @warn "scalar indexing into FieldVector (this can be very slow)" maxlog = 1 + getindex(fv, BlockArrays.findblockindex(axes(fv, 1), i)) +end +Base.@propagate_inbounds function Base.setindex!( + fv::Fields.FieldVector, + val, + i::Integer, +) + @warn "scalar indexing into FieldVector (this can be very slow)" maxlog = 1 + setindex!(fv, val, BlockArrays.findblockindex(axes(fv, 1), i)) +end + +# TODO: Move all of the following to ClimaCore.jl + +filtered_names(f::F, x) where {F} = filtered_names_at_name(f, x, @name()) +function filtered_names_at_name(f::F, x, name) where {F} + field = MatrixFields.get_field(x, name) + f(field) && return (name,) + internal_names = MatrixFields.top_level_names(field) + isempty(internal_names) && return () + tuples_of_names = MatrixFields.unrolled_map(internal_names) do internal_name + Base.@_inline_meta + child_name = MatrixFields.append_internal_name(name, internal_name) + filtered_names_at_name(f, x, child_name) + end + return MatrixFields.unrolled_flatten(tuples_of_names) +end +if hasfield(Method, :recursion_relation) + dont_limit = (args...) -> true + for m in methods(filtered_names_at_name) + m.recursion_relation = dont_limit + end +end + +scalar_field_names(fv) = + filtered_names(x -> x isa Fields.Field && eltype(x) == eltype(fv), fv) + +function scalar_level_iterator(field_vector) + Iterators.flatmap(scalar_field_names(field_vector)) do name + field = MatrixFields.get_field(field_vector, name) + if field isa Fields.SpectralElementField + (field,) + else + Iterators.map(1:Spaces.nlevels(axes(field))) do v + Fields.level(field, v - 1 + Operators.left_idx(axes(field))) + end + end + end +end + +function column_iterator(field_vector) + example_field = field_vector.c # TODO: Generalize this. + horz_space = Spaces.horizontal_space(axes(example_field)) + qs = 1:Quadratures.degrees_of_freedom(Spaces.quadrature_style(horz_space)) + hs = Spaces.eachslabindex(horz_space) + return if Fields.field_values(example_field) isa DataLayouts.VIFH + Iterators.map(Iterators.product(qs, hs)) do (i, h) + field_vector[Fields.ColumnIndex((i,), h)] + end + else + @assert Fields.field_values(example_field) isa DataLayouts.VIJFH + Iterators.map(Iterators.product(qs, qs, hs)) do (i, j, h) + field_vector[Fields.ColumnIndex((i, j), h)] + end + end +end diff --git a/src/prognostic_equations/implicit/exact_jacobian.jl b/src/prognostic_equations/implicit/exact_jacobian.jl new file mode 100644 index 0000000000..2e565313ad --- /dev/null +++ b/src/prognostic_equations/implicit/exact_jacobian.jl @@ -0,0 +1,196 @@ +@kwdef struct ExactJacobian <: JacobianAlgorithm + always_update_exact_jacobian::Bool = true + preserve_unfactorized_jacobian::Bool = false # only true for DebugJacobian +end + +function jacobian_cache(alg::ExactJacobian, Y, p) + FT = eltype(Y) + DA = ClimaComms.array_type(Y) + Y_columns = column_iterator(Y) + n_columns = length(Y_columns) + n_εs = length(first(Y_columns)) + dual_type = ForwardDiff.Dual{ExactJacobian, FT, 1} + + similar_with_dual_fields(named_tuple) = + Fields._values(similar(Fields.FieldVector(; named_tuple...), dual_type)) + + Y_dual = similar(Y, dual_type) + Y_dual_copy = similar(Y_dual) + Y_err_dual = similar(Y_dual) + # precomputed_dual = + # similar_with_dual_fields(precomputed_quantities(Y, p.atmos)) + implicit_precomputed_dual = + similar_with_dual_fields(implicit_precomputed_quantities(Y, p.atmos)) + scratch_dual = similar_with_dual_fields(p.scratch) + + column_partials_vector = DA{ForwardDiff.Partials{1, FT}}(undef, n_εs) + column_vectors = DA{FT}(undef, n_columns, n_εs) + column_matrices = DA{FT}(undef, n_columns, n_εs, n_εs) + column_factorized_matrices = + alg.preserve_unfactorized_jacobian ? copy(column_matrices) : + column_matrices + lu_cache = DA{FT}(undef, n_columns) + + return (; + Y_dual, + Y_dual_copy, + Y_err_dual, + implicit_precomputed_dual, # precomputed_dual, + scratch_dual, + column_partials_vector, + column_vectors, + column_matrices, + column_factorized_matrices, + lu_cache, + ) +end + +always_update_exact_jacobian(alg::ExactJacobian) = + alg.always_update_exact_jacobian + +# function factorize_exact_jacobian!(alg::ExactJacobian, cache, Y, p, dtγ, t) +# ... +# Y_dual .= Y +# zeros_tuple = ntuple(_ -> 0, n_εs) +# for (ε_index, Y_dual_scalar_level) in enumerate(Y_dual_scalar_levels) +# ε_partials = Base.setindex(zeros_tuple, 1, ε_index) +# ε = ForwardDiff.Dual{ExactJacobian}(0, ε_partials...) +# parent(Y_dual_scalar_level) .+= ε +# end # add a unique infinitesimal ε to the scalar values on each level of Y +# Y_dual_copy .= Y_dual # need a copy because Y_dual gets changed at boundary +# set_implicit_precomputed_quantities!(Y_dual, p_dual, t) # changes Y_dual +# implicit_tendency!(Y_err_dual, Y_dual, p_dual, t) # the tendency Jacobian +# @. Y_err_dual = dtγ * Y_err_dual - Y_dual_copy # the residual Jacobian +# for (col_index, Y_err_dual_column) in enumerate(Y_err_dual_columns) +# column_partials_vector = cache.column_partials_vectors[col_index] +# column_matrix = cache.column_matrices[col_index] +# column_factorized_matrix = cache.column_factorized_matrices[col_index] +# column_partials_vector .= ForwardDiff.partials.(Y_err_dual_column) +# column_matrix .= reinterpret(reshape, FT, column_partials_vector)' +# end # TODO: Parallelize this loop +# ... +# end +function factorize_exact_jacobian!(alg::ExactJacobian, cache, Y, p, dtγ, t) + (; preserve_unfactorized_jacobian) = alg + (; Y_dual, Y_dual_copy, Y_err_dual) = cache + # (; precomputed_dual, scratch_dual, column_partials_vector) = cache + (; implicit_precomputed_dual, scratch_dual, column_partials_vector) = cache + (; column_matrices, column_factorized_matrices, lu_cache) = cache + + set_precomputed_quantities!(Y, p, t) + p_dual_args = ntuple(Val(fieldcount(typeof(p)))) do cache_field_index + cache_field_name = fieldname(typeof(p), cache_field_index) + if cache_field_name == :precomputed + (; p.precomputed..., implicit_precomputed_dual...) # precomputed_dual + elseif cache_field_name == :scratch + scratch_dual + else + getfield(p, cache_field_index) + end + end + p_dual = AtmosCache(p_dual_args...) + + FT = eltype(Y) + Y_dual_scalar_levels = scalar_level_iterator(Y_dual) + Y_err_dual_columns = column_iterator(Y_err_dual) + ε = ForwardDiff.Dual{ExactJacobian}(0, 1) + + Y_dual .= Y + for (ε_index, Y_dual_scalar_level) in enumerate(Y_dual_scalar_levels) + parent(Y_dual_scalar_level) .+= ε # add ε to this level of Y_dual + Y_dual_copy .= Y_dual # make a copy because Y_dual gets changed at boundary + # set_precomputed_quantities!(Y_dual, p_dual, t) + set_implicit_precomputed_quantities!(Y_dual, p_dual, t) # changes Y_dual + implicit_tendency!(Y_err_dual, Y_dual, p_dual, t) # get tendency derivative + @. Y_err_dual = dtγ * Y_err_dual - Y_dual_copy # get residual derivative + parent(Y_dual_scalar_level) .-= ε # remove ε from this level of Y_dual + + for (col_index, Y_err_dual_column) in enumerate(Y_err_dual_columns) + column_partials_vector .= ForwardDiff.partials.(Y_err_dual_column) + view(column_matrices, col_index, :, ε_index) .= + reinterpret(reshape, FT, column_partials_vector) + end # TODO: Parallelize this loop. + end # TODO: Make this loop use batches of up to 32 values of ε at a time. + + first_column_matrix = view(cache.column_matrices, 1, :, :) + first_column_matrix_cpu = + ClimaComms.allowscalar(Array, ClimaComms.device(Y), first_column_matrix) + p.scratch.temp_matrix .= (first_column_matrix_cpu + I) / dtγ + + if preserve_unfactorized_jacobian + column_factorized_matrices .= column_matrices + end + parallel_lu_factorize!(column_factorized_matrices, lu_cache) +end + +approximate_jacobian!(::ExactJacobian, _, _, _, _, _) = nothing + +function invert_jacobian!(::ExactJacobian, cache, x, b) + (; column_vectors, column_factorized_matrices, lu_cache) = cache + for (col_index, b_column) in enumerate(column_iterator(b)) + view(column_vectors, col_index, :) .= b_column + end # TODO: Parallelize this loop. + parallel_lu_solve!(column_vectors, column_factorized_matrices, lu_cache) + for (col_index, x_column) in enumerate(column_iterator(x)) + x_column .= view(column_vectors, col_index, :) + end # TODO: Parallelize this loop. +end + +# Set the derivative of `sqrt(x)` to `iszero(x) ? zero(x) : inv(2 * sqrt(x))` in +# order to properly handle derivatives of `x * sqrt(x)`. Without this change, +# the derivative of `x * sqrt(x)` is `NaN` when `x` is zero. This method +# specializes on the tag `ExactJacobian` because not specializing on any tag +# overwrites the generic method for `Dual` in `ForwardDiff` and breaks +# precompilation, while specializing on the default tag `Nothing` causes the +# type piracy Aqua test to fail. +@inline function Base.sqrt(d::ForwardDiff.Dual{ExactJacobian}) + tag = Val{ExactJacobian}() + x = ForwardDiff.value(d) + partials = ForwardDiff.partials(d) + val = sqrt(x) + deriv = iszero(x) ? zero(x) : inv(2 * val) + return ForwardDiff.dual_definition_retval(tag, val, deriv, partials) +end + +function parallel_lu_factorize!(As, temporary_vector) + @assert ndims(As) == 3 + n = size(As, 2) + @assert size(As, 3) == n + @inbounds for k in 1:n + all(!isnan, view(As, :, k, k)) || error("LU error: NaN on diagonal") + all(!iszero, view(As, :, k, k)) || error("LU error: 0 on diagonal") + temporary_vector .= inv.(view(As, :, k, k)) + for i in (k + 1):n + view(As, :, i, k) .*= temporary_vector + end + for j in (k + 1):n + for i in (k + 1):n + view(As, :, i, j) .-= view(As, :, i, k) .* view(As, :, k, j) + end + end + end +end + +function parallel_lu_solve!(xs, As, temporary_vector) + @assert ndims(xs) == 2 && ndims(As) == 3 + n = size(As, 2) + @assert size(xs, 2) == n && size(As, 3) == n + @inbounds begin + for i in 2:n + temporary_vector .= zero(eltype(xs)) + for j in 1:(i - 1) + temporary_vector .+= view(As, :, i, j) .* view(xs, :, j) + end + view(xs, :, i) .-= temporary_vector + end + view(xs, :, n) ./= view(As, :, n, n) + for i in (n - 1):-1:1 + temporary_vector .= zero(eltype(xs)) + for j in (i + 1):n + temporary_vector .+= view(As, :, i, j) .* view(xs, :, j) + end + view(xs, :, i) .-= temporary_vector + view(xs, :, i) ./= view(As, :, i, i) + end + end +end diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl index d8f721e967..c818aeb0fd 100644 --- a/src/prognostic_equations/implicit/implicit_solver.jl +++ b/src/prognostic_equations/implicit/implicit_solver.jl @@ -1,28 +1,25 @@ -import LinearAlgebra: I, Adjoint, ldiv! +import BlockArrays, ForwardDiff +import LinearAlgebra: I, Adjoint, LU, lu, lu!, ldiv! import ClimaCore.MatrixFields: @name -using ClimaCore.MatrixFields - -abstract type DerivativeFlag end -struct UseDerivative <: DerivativeFlag end -struct IgnoreDerivative <: DerivativeFlag end -use_derivative(::UseDerivative) = true -use_derivative(::IgnoreDerivative) = false """ - ImplicitEquationJacobian( - Y, atmos; - approximate_solve_iters, diffusion_flag, topography_flag, sgs_advection_flag, transform_flag - ) + JacobianAlgorithm -A wrapper for the matrix ``∂E/∂Y``, where ``E(Y)`` is the "error" of the -implicit step with the state ``Y``. +A description of how to compute the matrix ``∂R/∂Y``, where ``R(Y)`` denotes the +residual of an implicit step with the state ``Y``. Concrete implementations of +this abstract type should define 4 methods: + - `jacobian_cache(alg::JacobianAlgorithm, Y, p)` + - `always_update_exact_jacobian(alg::JacobianAlgorithm)` + - `factorize_exact_jacobian!(alg::JacobianAlgorithm, cache, Y, p, dtγ, t)` + - `approximate_jacobian!(alg::JacobianAlgorithm, cache, Y, p, dtγ, t)` + - `invert_jacobian!(alg::JacobianAlgorithm, cache, x, b)` # Background When we use an implicit or split implicit-explicit (IMEX) timestepping scheme, -we end up with a nonlinear equation of the form ``E(Y) = 0``, where +we end up with a nonlinear equation of the form ``R(Y) = 0``, where ```math - E(Y) = Y_{imp}(Y) - Y = \\hat{Y} + Δt * T_{imp}(Y) - Y. + R(Y) = Y_{imp}(Y) - Y = \\hat{Y} + Δt * T_{imp}(Y) - Y. ``` In this expression, ``Y_{imp}(Y)`` denotes the state at some time ``t + Δt``. This can be expressed as the sum of ``\\hat{Y}``, the contribution from the @@ -38,827 +35,90 @@ divided into several sub-steps or "stages", where the duration of stage ``i`` is ``Δt * γ_i`` for some constant ``γ_i`` between 0 and 1. In order to solve this equation using Newton's method, we must specify the -derivative ``∂E/∂Y``. Since ``\\hat{Y}`` does not depend on ``Y`` (it is only a +derivative ``∂R/∂Y``. Since ``\\hat{Y}`` does not depend on ``Y`` (it is only a function of the state at or before time ``t``), this derivative is ```math - E'(Y) = Δt * T_{imp}'(Y) - I. + R'(Y) = Δt * T_{imp}'(Y) - I. ``` -In addition, we must specify how to divide ``E(Y)`` by this derivative, i.e., +In addition, we must specify how to divide ``R(Y)`` by this derivative, i.e., how to solve the linear equation ```math - E'(Y) * ΔY = E(Y). + R'(Y) * ΔY = R(Y). ``` Note: This equation comes from assuming that there is some ``ΔY`` such that -``E(Y - ΔY) = 0`` and making the first-order approximation +``R(Y - ΔY) = 0`` and making the first-order approximation ```math - E(Y - ΔY) \\approx E(Y) - E'(Y) * ΔY. + R(Y - ΔY) \\approx R(Y) - R'(Y) * ΔY. ``` After initializing ``Y`` to ``Y[0] = \\hat{Y}``, Newton's method executes the following steps: -- Compute the derivative ``E'(Y[0])``. -- Compute the implicit tendency ``T_{imp}(Y[0])`` and use it to get ``E(Y[0])``. -- Solve the linear equation ``E'(Y[0]) * ΔY[0] = E(Y[0])`` for ``ΔY[0]``. +- Compute the derivative ``R'(Y[0])``. +- Compute the implicit tendency ``T_{imp}(Y[0])`` and use it to get ``R(Y[0])``. +- Solve the linear equation ``R'(Y[0]) * ΔY[0] = R(Y[0])`` for ``ΔY[0]``. - Update ``Y`` to ``Y[1] = Y[0] - ΔY[0]``. If the number of Newton iterations is limited to 1, this new value of ``Y`` is taken to be the solution of the implicit equation. Otherwise, this sequence of steps is repeated, i.e., ``ΔY[1]`` is computed and used to update ``Y`` to ``Y[2] = Y[1] - ΔY[1]``, then ``ΔY[2]`` is computed and used to update ``Y`` to ``Y[3] = Y[2] - ΔY[2]``, and so on. The iterative process is terminated either -when the error ``E(Y)`` is sufficiently close to 0 (according to the convergence -condition passed to Newton's method), or when the maximum number of iterations -is reached. - -# Arguments - -- `Y::FieldVector`: the state of the simulation -- `atmos::AtmosModel`: the model configuration -- `approximate_solve_iters::Int`: number of iterations to take for the - approximate linear solve required when `diffusion_flag = UseDerivative()` -- `diffusion_flag::DerivativeFlag`: whether the derivative of the - diffusion tendency with respect to the quantities being diffused should be - computed or approximated as 0; must be either `UseDerivative()` or - `IgnoreDerivative()` instead of a `Bool` to ensure type-stability -- `topography_flag::DerivativeFlag`: whether the derivative of vertical - contravariant velocity with respect to horizontal covariant velocity should be - computed or approximated as 0; must be either `UseDerivative()` or - `IgnoreDerivative()` instead of a `Bool` to ensure type-stability -- `sgs_advection_flag::DerivativeFlag`: whether the derivative of the - subgrid-scale advection tendency with respect to the updraft quantities should be - computed or approximated as 0; must be either `UseDerivative()` or - `IgnoreDerivative()` instead of a `Bool` to ensure type-stability -- `transform_flag::Bool`: whether the error should be transformed from ``E(Y)`` - to ``E(Y)/Δt``, which is required for non-Rosenbrock timestepping schemes from - OrdinaryDiffEq.jl +when the residual ``R(Y)`` is sufficiently close to 0 (according to the +convergence condition passed to Newton's method), or when the maximum number of +iterations is reached. """ -struct ImplicitEquationJacobian{ - M <: MatrixFields.FieldMatrix, - S <: MatrixFields.FieldMatrixSolver, - F1 <: DerivativeFlag, - F2 <: DerivativeFlag, - F3 <: DerivativeFlag, - T <: Fields.FieldVector, - R <: Base.RefValue, -} - # stores the matrix E'(Y) = Δt * T_imp'(Y) - I - matrix::M - - # solves the linear equation E'(Y) * ΔY = E(Y) for ΔY - solver::S +abstract type JacobianAlgorithm end - # flags that determine how E'(Y) is approximated - diffusion_flag::F1 - topography_flag::F2 - sgs_advection_flag::F3 - - # required by Krylov.jl to evaluate ldiv! with AbstractVector inputs - temp_b::T - temp_x::T - - # required by OrdinaryDiffEq.jl to run non-Rosenbrock timestepping schemes - transform_flag::Bool - dtγ_ref::R +struct ImplicitEquationJacobian{A <: JacobianAlgorithm, C} + alg::A + cache::C end - -function ImplicitEquationJacobian( - Y, - atmos; - approximate_solve_iters = 1, - diffusion_flag = atmos.diff_mode == Implicit() ? UseDerivative() : - IgnoreDerivative(), - topography_flag = has_topography(axes(Y.c)) ? UseDerivative() : - IgnoreDerivative(), - sgs_advection_flag = atmos.sgs_adv_mode == Implicit() ? UseDerivative() : - IgnoreDerivative(), - transform_flag = false, -) - FT = Spaces.undertype(axes(Y.c)) - CTh = CTh_vector_type(axes(Y.c)) - - DiagonalRow = DiagonalMatrixRow{FT} - TridiagonalRow = TridiagonalMatrixRow{FT} - BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}} - TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}} - BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}} - BidiagonalRow_C3xACTh = - BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CTh{FT})')} - DiagonalRow_C3xACT3 = - DiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} - TridiagonalRow_C3xACT3 = - TridiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} - - is_in_Y(name) = MatrixFields.has_field(Y, name) - - ρq_tot_if_available = is_in_Y(@name(c.ρq_tot)) ? (@name(c.ρq_tot),) : () - ρatke_if_available = - is_in_Y(@name(c.sgs⁰.ρatke)) ? (@name(c.sgs⁰.ρatke),) : () - sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : () - - tracer_names = ( - @name(c.ρq_tot), - @name(c.ρq_liq), - @name(c.ρq_ice), - @name(c.ρq_rai), - @name(c.ρq_sno), - ) - available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names) - - # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0, - # which means that multiplying inv(-1) by a Float32 will yield a Float64. - identity_blocks = MatrixFields.unrolled_map( - name -> (name, name) => FT(-1) * I, - (@name(c.ρ), sfc_if_available...), +function ImplicitEquationJacobian(alg, Y, p) + cache = (; + jacobian_cache(alg, Y, p)..., + x_krylov = similar(Y), + b_krylov = similar(Y), ) + return ImplicitEquationJacobian(alg, cache) +end - active_scalar_names = (@name(c.ρ), @name(c.ρe_tot), ρq_tot_if_available...) - advection_blocks = ( - ( - use_derivative(topography_flag) ? - MatrixFields.unrolled_map( - name -> - (name, @name(c.uₕ)) => - similar(Y.c, TridiagonalRow_ACTh), - active_scalar_names, - ) : () - )..., - MatrixFields.unrolled_map( - name -> (name, @name(f.u₃)) => similar(Y.c, BidiagonalRow_ACT3), - active_scalar_names, - )..., - MatrixFields.unrolled_map( - name -> (@name(f.u₃), name) => similar(Y.f, BidiagonalRow_C3), - active_scalar_names, - )..., - (@name(f.u₃), @name(c.uₕ)) => similar(Y.f, BidiagonalRow_C3xACTh), - (@name(f.u₃), @name(f.u₃)) => similar(Y.f, TridiagonalRow_C3xACT3), - ) - - diffused_scalar_names = - (@name(c.ρe_tot), available_tracer_names..., ρatke_if_available...) - diffusion_blocks = if use_derivative(diffusion_flag) - ( - MatrixFields.unrolled_map( - name -> (name, @name(c.ρ)) => similar(Y.c, TridiagonalRow), - diffused_scalar_names, - )..., - MatrixFields.unrolled_map( - name -> (name, name) => similar(Y.c, TridiagonalRow), - diffused_scalar_names, - )..., - ( - is_in_Y(@name(c.ρq_tot)) ? - ( - (@name(c.ρe_tot), @name(c.ρq_tot)) => - similar(Y.c, TridiagonalRow), - ) : () - )..., - (@name(c.uₕ), @name(c.uₕ)) => - !isnothing(atmos.turbconv_model) || - diffuse_momentum(atmos.vert_diff) ? - similar(Y.c, TridiagonalRow) : FT(-1) * I, - ) - else - MatrixFields.unrolled_map( - name -> (name, name) => FT(-1) * I, - (diffused_scalar_names..., @name(c.uₕ)), - ) - end - - sgs_advection_blocks = if atmos.turbconv_model isa PrognosticEDMFX - @assert n_prognostic_mass_flux_subdomains(atmos.turbconv_model) == 1 - sgs_scalar_names = ( - @name(c.sgsʲs.:(1).q_tot), - @name(c.sgsʲs.:(1).mse), - @name(c.sgsʲs.:(1).ρa) - ) - if use_derivative(sgs_advection_flag) - ( - MatrixFields.unrolled_map( - name -> (name, name) => similar(Y.c, TridiagonalRow), - sgs_scalar_names, - )..., - (@name(c.sgsʲs.:(1).mse), @name(c.ρ)) => - similar(Y.c, DiagonalRow), - (@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)) => - similar(Y.c, DiagonalRow), - (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)) => - similar(Y.c, TridiagonalRow), - (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)) => - similar(Y.c, TridiagonalRow), - (@name(f.sgsʲs.:(1).u₃), @name(c.ρ)) => - similar(Y.f, BidiagonalRow_C3), - (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)) => - similar(Y.f, BidiagonalRow_C3), - (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)) => - similar(Y.f, BidiagonalRow_C3), - (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => - similar(Y.f, TridiagonalRow_C3xACT3), - ) - else - ( - MatrixFields.unrolled_map( - name -> (name, name) => FT(-1) * I, - sgs_scalar_names, - )..., - (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => - !isnothing(atmos.rayleigh_sponge) ? - similar(Y.f, DiagonalRow_C3xACT3) : FT(-1) * I, - ) - end - else - () - end +# We only use A, but ClimaTimeSteppers.jl currently requires us to pass +# jac_prototype and then calls similar(jac_prototype) to obtain A. Since +# this was only done for compatibility with OrdinaryDiffEq.jl, which we are no +# longer using, we should fix this in a new version of ClimaTimeSteppers.jl. +Base.similar(A::ImplicitEquationJacobian) = A - matrix = MatrixFields.FieldMatrix( - identity_blocks..., - sgs_advection_blocks..., - advection_blocks..., - diffusion_blocks..., - ) +# This is called from a callback when always_update_exact_jacobian is false. +NVTX.@annotate function update_exact_jacobian!(A, Y, p, dtγ, t) + FT = eltype(Y) + factorize_exact_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t) +end - sgs_names_if_available = if atmos.turbconv_model isa PrognosticEDMFX - ( - @name(c.sgsʲs.:(1).q_tot), - @name(c.sgsʲs.:(1).mse), - @name(c.sgsʲs.:(1).ρa), - @name(f.sgsʲs.:(1).u₃), - ) - else - () +# This is passed to ClimaTimeSteppers.jl and called on each Newton iteration. +NVTX.@annotate function update_jacobian!(A, Y, p, dtγ, t) + FT = eltype(Y) + if always_update_exact_jacobian(A.alg) + factorize_exact_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t) end - names₁_group₁ = (@name(c.ρ), sfc_if_available...) - names₁_group₂ = (available_tracer_names..., ρatke_if_available...) - names₁_group₃ = (@name(c.ρe_tot),) - names₁ = ( - names₁_group₁..., - sgs_names_if_available..., - names₁_group₂..., - names₁_group₃..., - ) - - alg₂ = MatrixFields.BlockLowerTriangularSolve(@name(c.uₕ)) - alg = - if use_derivative(diffusion_flag) || use_derivative(sgs_advection_flag) - alg₁_subalg₂ = - if atmos.turbconv_model isa PrognosticEDMFX && - use_derivative(sgs_advection_flag) - diff_subalg = - use_derivative(diffusion_flag) ? - (; - alg₂ = MatrixFields.BlockLowerTriangularSolve( - names₁_group₂..., - ) - ) : (;) - (; - alg₂ = MatrixFields.BlockLowerTriangularSolve( - @name(c.sgsʲs.:(1).q_tot); - alg₂ = MatrixFields.BlockLowerTriangularSolve( - @name(c.sgsʲs.:(1).mse); - alg₂ = MatrixFields.BlockLowerTriangularSolve( - @name(c.sgsʲs.:(1).ρa), - @name(f.sgsʲs.:(1).u₃); - diff_subalg..., - ), - ), - ) - ) - else - is_in_Y(@name(c.ρq_tot)) ? - (; - alg₂ = MatrixFields.BlockLowerTriangularSolve( - names₁_group₂..., - ) - ) : (;) - end - alg₁ = MatrixFields.BlockLowerTriangularSolve( - names₁_group₁...; - alg₁_subalg₂..., - ) - MatrixFields.ApproximateBlockArrowheadIterativeSolve( - names₁...; - alg₁, - alg₂, - P_alg₁ = MatrixFields.MainDiagonalPreconditioner(), - n_iters = approximate_solve_iters, - ) - else - MatrixFields.BlockArrowheadSolve(names₁...; alg₂) - end - - return ImplicitEquationJacobian( - matrix, - MatrixFields.FieldMatrixSolver(alg, matrix, Y), - diffusion_flag, - topography_flag, - sgs_advection_flag, - similar(Y), - similar(Y), - transform_flag, - Ref{FT}(), - ) + approximate_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t) end -# We only use A, but ClimaTimeSteppers.jl require us to -# pass jac_prototype and then call similar(jac_prototype) to -# obtain A. This is a workaround to avoid unnecessary allocations. -Base.similar(A::ImplicitEquationJacobian) = A - -# This method specifies how to solve the equation E'(Y) * ΔY = E(Y) for ΔY. -NVTX.@annotate function ldiv!( +# This is called by ClimaTimeSteppers.jl on each Newton iteration. +NVTX.@annotate ldiv!( x::Fields.FieldVector, A::ImplicitEquationJacobian, b::Fields.FieldVector, -) - MatrixFields.field_matrix_solve!(A.solver, x, A.matrix, b) - if A.transform_flag - @. x *= -A.dtγ_ref[] - end -end +) = invert_jacobian!(A.alg, A.cache, x, b) -# This method for ldiv! is called by Krylov.jl from inside ClimaTimeSteppers.jl. -# See https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/605 for a -# related issue that requires the same workaround. -NVTX.@annotate function ldiv!( +# This is called by Krylov.jl from inside ClimaTimeSteppers.jl. See +# https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/605 for a related +# issue that requires the same workaround. +function ldiv!( x::AbstractVector, A::ImplicitEquationJacobian, b::AbstractVector, ) - A.temp_b .= b - ldiv!(A.temp_x, A, A.temp_b) - x .= A.temp_x -end - -# This function is used by DiffEqBase.jl instead of ldiv!. -linsolve!(::Type{Val{:init}}, f, u0; kwargs...) = _linsolve! -_linsolve!(x, A, b, update_matrix = false; kwargs...) = ldiv!(x, A, b) - -# This method specifies how to compute E'(Y), which is referred to as "Wfact" in -# DiffEqBase.jl. -NVTX.@annotate function Wfact!(A, Y, p, dtγ, t) - # Remove unnecessary values from p to avoid allocations in bycolumn. - p′ = (; - p.precomputed.ᶜspecific, - p.precomputed.ᶜK, - p.precomputed.ᶜts, - p.precomputed.ᶜp, - ( - p.atmos.precip_model isa Microphysics1Moment ? - (; p.precomputed.ᶜwᵣ, p.precomputed.ᶜwₛ) : (;) - )..., - p.precomputed.ᶜh_tot, - ( - use_derivative(A.diffusion_flag) ? - (; p.precomputed.ᶜK_u, p.precomputed.ᶜK_h) : (;) - )..., - ( - use_derivative(A.diffusion_flag) && - p.atmos.turbconv_model isa AbstractEDMF ? - (; p.precomputed.ᶜtke⁰, p.precomputed.ᶜmixing_length) : (;) - )..., - ( - use_derivative(A.diffusion_flag) && - p.atmos.turbconv_model isa PrognosticEDMFX ? - (; p.precomputed.ᶜρa⁰) : (;) - )..., - ( - use_derivative(A.sgs_advection_flag) && - p.atmos.turbconv_model isa PrognosticEDMFX ? - (; - p.core.ᶜgradᵥ_ᶠΦ, - p.precomputed.ᶜρʲs, - p.precomputed.ᶠu³ʲs, - p.precomputed.ᶜtsʲs, - p.precomputed.bdmr_l, - p.precomputed.bdmr_r, - p.precomputed.bdmr, - ) : (;) - )..., - p.core.ᶜΦ, - p.core.ᶠgradᵥ_ᶜΦ, - p.scratch.ᶜtemp_scalar, - p.scratch.ᶜtemp_C3, - p.scratch.ᶠtemp_CT3, - p.scratch.∂ᶜK_∂ᶜuₕ, - p.scratch.∂ᶜK_∂ᶠu₃, - p.scratch.ᶠp_grad_matrix, - p.scratch.ᶜadvection_matrix, - p.scratch.ᶜdiffusion_h_matrix, - p.scratch.ᶜdiffusion_u_matrix, - p.scratch.ᶠbidiagonal_matrix_ct3, - p.scratch.ᶠbidiagonal_matrix_ct3_2, - p.scratch.ᶠtridiagonal_matrix_c3, - p.dt, - p.params, - p.atmos, - ( - p.atmos.rayleigh_sponge isa RayleighSponge ? - (; p.rayleigh_sponge.ᶠβ_rayleigh_w) : (;) - )..., - ) - - # Convert dtγ from a Float64 to an FT. - FT = Spaces.undertype(axes(Y.c)) - dtγ′ = FT(dtγ) - - A.dtγ_ref[] = dtγ′ - update_implicit_equation_jacobian!(A, Y, p′, dtγ′) -end - -function update_implicit_equation_jacobian!(A, Y, p, dtγ) - (; matrix, diffusion_flag, sgs_advection_flag, topography_flag) = A - (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶜΦ, ᶠgradᵥ_ᶜΦ) = p - (; - ᶜtemp_C3, - ∂ᶜK_∂ᶜuₕ, - ∂ᶜK_∂ᶠu₃, - ᶠp_grad_matrix, - ᶜadvection_matrix, - ᶠbidiagonal_matrix_ct3, - ᶠbidiagonal_matrix_ct3_2, - ᶠtridiagonal_matrix_c3, - ) = p - (; ᶜdiffusion_h_matrix, ᶜdiffusion_u_matrix, params) = p - (; edmfx_upwinding) = p.atmos.numerics - - FT = Spaces.undertype(axes(Y.c)) - CTh = CTh_vector_type(axes(Y.c)) - one_C3xACT3 = C3(FT(1)) * CT3(FT(1))' - - cv_d = FT(CAP.cv_d(params)) - Δcv_v = FT(CAP.cv_v(params)) - cv_d - T_tri = FT(CAP.T_triple(params)) - e_int_v0 = T_tri * Δcv_v - FT(CAP.e_int_v0(params)) - thermo_params = CAP.thermodynamics_params(params) - - ᶜρ = Y.c.ρ - ᶜuₕ = Y.c.uₕ - ᶠu₃ = Y.f.u₃ - ᶜJ = Fields.local_geometry_field(Y.c).J - ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ - ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ - - ᶜkappa_m = p.ᶜtemp_scalar - @. ᶜkappa_m = - TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts) - - if use_derivative(topography_flag) - @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow( - adjoint(CTh(ᶜuₕ)) + adjoint(ᶜinterp(ᶠu₃)) * g³ʰ(ᶜgⁱʲ), - ) - else - @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ))) - end - @. ∂ᶜK_∂ᶠu₃ = - ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) + - DiagonalMatrixRow(adjoint(CT3(ᶜuₕ))) ⋅ ᶜinterp_matrix() - - @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() - - @. ᶜadvection_matrix = - -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) - - if use_derivative(topography_flag) - ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)] - @. ∂ᶜρ_err_∂ᶜuₕ = - dtγ * ᶜadvection_matrix ⋅ ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ - DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ)) - end - ∂ᶜρ_err_∂ᶠu₃ = matrix[@name(c.ρ), @name(f.u₃)] - @. ∂ᶜρ_err_∂ᶠu₃ = dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ)) - - tracer_info = ( - (@name(c.ρe_tot), @name(ᶜh_tot)), - (@name(c.ρq_tot), @name(ᶜspecific.q_tot)), - ) - MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name) - MatrixFields.has_field(Y, ρχ_name) || return - ᶜχ = MatrixFields.get_field(p, χ_name) - if use_derivative(topography_flag) - ∂ᶜρχ_err_∂ᶜuₕ = matrix[ρχ_name, @name(c.uₕ)] - end - ∂ᶜρχ_err_∂ᶠu₃ = matrix[ρχ_name, @name(f.u₃)] - use_derivative(topography_flag) && @. ∂ᶜρχ_err_∂ᶜuₕ = - dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅ - ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ)) - @. ∂ᶜρχ_err_∂ᶠu₃ = - dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ) * g³³(ᶠgⁱʲ)) - end - - ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)] - ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)] - @. ∂ᶠu₃_err_∂ᶜρ = - dtγ * ( - ᶠp_grad_matrix ⋅ - DiagonalMatrixRow(ᶜkappa_m * (T_tri * cv_d - ᶜK - ᶜΦ)) + - DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅ - ᶠinterp_matrix() - ) - @. ∂ᶠu₃_err_∂ᶜρe_tot = dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m) - if MatrixFields.has_field(Y, @name(c.ρq_tot)) - ∂ᶠu₃_err_∂ᶜρq_tot = matrix[@name(f.u₃), @name(c.ρq_tot)] - @. ∂ᶠu₃_err_∂ᶜρq_tot = - dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m * e_int_v0) - end - - ∂ᶠu₃_err_∂ᶜuₕ = matrix[@name(f.u₃), @name(c.uₕ)] - ∂ᶠu₃_err_∂ᶠu₃ = matrix[@name(f.u₃), @name(f.u₃)] - I_u₃ = DiagonalMatrixRow(one_C3xACT3) - @. ∂ᶠu₃_err_∂ᶜuₕ = - dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ ∂ᶜK_∂ᶜuₕ - if p.atmos.rayleigh_sponge isa RayleighSponge - @. ∂ᶠu₃_err_∂ᶠu₃ = - dtγ * ( - ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ - ∂ᶜK_∂ᶠu₃ + DiagonalMatrixRow(-p.ᶠβ_rayleigh_w * (one_C3xACT3,)) - ) - (I_u₃,) - else - @. ∂ᶠu₃_err_∂ᶠu₃ = - dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ - ∂ᶜK_∂ᶠu₃ - (I_u₃,) - end - - if use_derivative(diffusion_flag) - (; ᶜK_h, ᶜK_u) = p - @. ᶜdiffusion_h_matrix = - ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅ - ᶠgradᵥ_matrix() - if ( - MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) || - !isnothing(p.atmos.turbconv_model) || - diffuse_momentum(p.atmos.vert_diff) - ) - @. ᶜdiffusion_u_matrix = - ᶜadvdivᵥ_matrix() ⋅ - DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)) ⋅ ᶠgradᵥ_matrix() - end - - ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name(c.ρe_tot), @name(c.ρ)] - ∂ᶜρe_tot_err_∂ᶜρe_tot = matrix[@name(c.ρe_tot), @name(c.ρe_tot)] - @. ∂ᶜρe_tot_err_∂ᶜρ = - dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow( - ( - -(1 + ᶜkappa_m) * ᶜspecific.e_tot - - ᶜkappa_m * e_int_v0 * ᶜspecific.q_tot - ) / ᶜρ, - ) - @. ∂ᶜρe_tot_err_∂ᶜρe_tot = - dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow((1 + ᶜkappa_m) / ᶜρ) - - (I,) - if MatrixFields.has_field(Y, @name(c.ρq_tot)) - ∂ᶜρe_tot_err_∂ᶜρq_tot = matrix[@name(c.ρe_tot), @name(c.ρq_tot)] - @. ∂ᶜρe_tot_err_∂ᶜρq_tot = - dtγ * ᶜdiffusion_h_matrix ⋅ - DiagonalMatrixRow(ᶜkappa_m * e_int_v0 / ᶜρ) - end - - tracer_info = ( - (@name(c.ρq_tot), @name(q_tot)), - (@name(c.ρq_liq), @name(q_liq)), - (@name(c.ρq_ice), @name(q_ice)), - (@name(c.ρq_rai), @name(q_rai)), - (@name(c.ρq_sno), @name(q_sno)), - ) - MatrixFields.unrolled_foreach(tracer_info) do (ρq_name, q_name) - MatrixFields.has_field(Y, ρq_name) || return - ᶜq = MatrixFields.get_field(ᶜspecific, q_name) - ∂ᶜρq_err_∂ᶜρ = matrix[ρq_name, @name(c.ρ)] - ∂ᶜρq_err_∂ᶜρq = matrix[ρq_name, ρq_name] - @. ∂ᶜρq_err_∂ᶜρ = - dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ) - @. ∂ᶜρq_err_∂ᶜρq = - dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ) - (I,) - end - - if MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) - turbconv_params = CAP.turbconv_params(params) - c_d = CAP.tke_diss_coeff(turbconv_params) - (; ᶜtke⁰, ᶜmixing_length, dt) = p - ᶜρa⁰ = p.atmos.turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : ᶜρ - ᶜρatke⁰ = Y.c.sgs⁰.ρatke - - @inline dissipation_rate(tke⁰, mixing_length) = - tke⁰ >= 0 ? c_d * sqrt(tke⁰) / max(mixing_length, 1) : 1 / dt - @inline ∂dissipation_rate_∂tke⁰(tke⁰, mixing_length) = - tke⁰ > 0 ? c_d / (2 * max(mixing_length, 1) * sqrt(tke⁰)) : - typeof(tke⁰)(0) - - ᶜdissipation_matrix_diagonal = p.ᶜtemp_scalar - @. ᶜdissipation_matrix_diagonal = - ᶜρatke⁰ * ∂dissipation_rate_∂tke⁰(ᶜtke⁰, ᶜmixing_length) - - ∂ᶜρatke⁰_err_∂ᶜρ = matrix[@name(c.sgs⁰.ρatke), @name(c.ρ)] - ∂ᶜρatke⁰_err_∂ᶜρatke⁰ = - matrix[@name(c.sgs⁰.ρatke), @name(c.sgs⁰.ρatke)] - @. ∂ᶜρatke⁰_err_∂ᶜρ = - dtγ * ( - ᶜdiffusion_u_matrix - - DiagonalMatrixRow(ᶜdissipation_matrix_diagonal) - ) ⋅ DiagonalMatrixRow(-(ᶜtke⁰) / ᶜρa⁰) - @. ∂ᶜρatke⁰_err_∂ᶜρatke⁰ = - dtγ * ( - ( - ᶜdiffusion_u_matrix - - DiagonalMatrixRow(ᶜdissipation_matrix_diagonal) - ) ⋅ DiagonalMatrixRow(1 / ᶜρa⁰) - - DiagonalMatrixRow(dissipation_rate(ᶜtke⁰, ᶜmixing_length)) - ) - (I,) - end - - if ( - !isnothing(p.atmos.turbconv_model) || - diffuse_momentum(p.atmos.vert_diff) - ) - ∂ᶜuₕ_err_∂ᶜuₕ = matrix[@name(c.uₕ), @name(c.uₕ)] - @. ∂ᶜuₕ_err_∂ᶜuₕ = - dtγ * DiagonalMatrixRow(1 / ᶜρ) ⋅ ᶜdiffusion_u_matrix - (I,) - end - - ᶠlg = Fields.local_geometry_field(Y.f) - precip_info = - ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ))) - MatrixFields.unrolled_foreach(precip_info) do (ρqₚ_name, wₚ_name) - MatrixFields.has_field(Y, ρqₚ_name) || return - ∂ᶜρqₚ_err_∂ᶜρqₚ = matrix[ρqₚ_name, ρqₚ_name] - ᶜwₚ = MatrixFields.get_field(p, wₚ_name) - ᶠtmp = p.ᶠtemp_CT3 - @. ᶠtmp = CT3(unit_basis_vector_data(CT3, ᶠlg)) * ᶠwinterp(ᶜJ, ᶜρ) - @. ∂ᶜρqₚ_err_∂ᶜρqₚ += - dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅ - ᶠright_bias_matrix() ⋅ DiagonalMatrixRow(-(ᶜwₚ) / ᶜρ) - end - end - - if p.atmos.turbconv_model isa PrognosticEDMFX - if use_derivative(sgs_advection_flag) - (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p - (; bdmr_l, bdmr_r, bdmr) = p - is_third_order = edmfx_upwinding == Val(:third_order) - ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1 - ᶠset_upwind_bcs = Operators.SetBoundaryOperator(; - top = Operators.SetValue(zero(CT3{FT})), - bottom = Operators.SetValue(zero(CT3{FT})), - ) # Need to wrap ᶠupwind in this for well-defined boundaries. - UpwindMatrixRowType = - is_third_order ? QuaddiagonalMatrixRow : BidiagonalMatrixRow - ᶠupwind_matrix = is_third_order ? ᶠupwind3_matrix : ᶠupwind1_matrix - ᶠset_upwind_matrix_bcs = Operators.SetBoundaryOperator(; - top = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), - bottom = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), - ) # Need to wrap ᶠupwind_matrix in this for well-defined boundaries. - ᶜkappa_mʲ = p.ᶜtemp_scalar - @. ᶜkappa_mʲ = - TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)) / - TD.cv_m(thermo_params, ᶜtsʲs.:(1)) - ∂ᶜq_totʲ_err_∂ᶜq_totʲ = - matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)] - @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ = - dtγ * ( - DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) - - ᶜadvdivᵥ_matrix() ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) - ) - (I,) - - ∂ᶜmseʲ_err_∂ᶜq_totʲ = - matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)] - @. ∂ᶜmseʲ_err_∂ᶜq_totʲ = - dtγ * ( - -DiagonalMatrixRow( - adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * - ᶜgradᵥ_ᶠΦ * - Y.c.ρ * - ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * e_int_v0, - ) - ) - ∂ᶜmseʲ_err_∂ᶜρ = matrix[@name(c.sgsʲs.:(1).mse), @name(c.ρ)] - @. ∂ᶜmseʲ_err_∂ᶜρ = - dtγ * ( - -DiagonalMatrixRow( - adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * ᶜgradᵥ_ᶠΦ / ᶜρʲs.:(1), - ) - ) - ∂ᶜmseʲ_err_∂ᶜmseʲ = - matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).mse)] - @. ∂ᶜmseʲ_err_∂ᶜmseʲ = - dtγ * ( - DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) - - ᶜadvdivᵥ_matrix() ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) - - DiagonalMatrixRow( - adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * - ᶜgradᵥ_ᶠΦ * - Y.c.ρ * - ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp), - ) - ) - (I,) - - ∂ᶜρaʲ_err_∂ᶜq_totʲ = - matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)] - @. ᶠbidiagonal_matrix_ct3 = - DiagonalMatrixRow( - ᶠset_upwind_bcs( - ᶠupwind( - ᶠu³ʲs.:(1), - draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)), - ), - ), - ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) * - e_int_v0, - ) - @. ᶠbidiagonal_matrix_ct3_2 = - DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ - DiagonalMatrixRow( - Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * - e_int_v0, - ) - - @. ∂ᶜρaʲ_err_∂ᶜq_totʲ = - dtγ * ᶜadvdivᵥ_matrix() ⋅ - (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2) - ∂ᶜρaʲ_err_∂ᶜmseʲ = - matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)] - @. ᶠbidiagonal_matrix_ct3 = - DiagonalMatrixRow( - ᶠset_upwind_bcs( - ᶠupwind( - ᶠu³ʲs.:(1), - draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)), - ), - ), - ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp), - ) - @. ᶠbidiagonal_matrix_ct3_2 = - DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ - DiagonalMatrixRow( - Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp), - ) - @. ∂ᶜρaʲ_err_∂ᶜmseʲ = - dtγ * ᶜadvdivᵥ_matrix() ⋅ - (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2) - ∂ᶜρaʲ_err_∂ᶜρaʲ = - matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)] - @. ᶜadvection_matrix = - -(ᶜadvdivᵥ_matrix()) ⋅ - DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) - @. ∂ᶜρaʲ_err_∂ᶜρaʲ = - dtγ * ᶜadvection_matrix ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅ - DiagonalMatrixRow(1 / ᶜρʲs.:(1)) - (I,) - - ∂ᶠu₃ʲ_err_∂ᶜρ = matrix[@name(f.sgsʲs.:(1).u₃), @name(c.ρ)] - @. ∂ᶠu₃ʲ_err_∂ᶜρ = - dtγ * DiagonalMatrixRow(ᶠgradᵥ_ᶜΦ / ᶠinterp(ᶜρʲs.:(1))) ⋅ - ᶠinterp_matrix() - ∂ᶠu₃ʲ_err_∂ᶜq_totʲ = - matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)] - @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ = - dtγ * DiagonalMatrixRow( - ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2, - ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) * - e_int_v0, - ) - ∂ᶠu₃ʲ_err_∂ᶜmseʲ = - matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)] - @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ = - dtγ * DiagonalMatrixRow( - ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2, - ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp), - ) - ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = - matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] - ᶜu₃ʲ = ᶜtemp_C3 - @. ᶜu₃ʲ = ᶜinterp(Y.f.sgsʲs.:(1).u₃) - - @. bdmr_l = convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix()) - @. bdmr_r = convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix()) - @. bdmr = ifelse(ᶜu₃ʲ.components.data.:1 > 0, bdmr_l, bdmr_r) - - @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋅ bdmr - if p.atmos.rayleigh_sponge isa RayleighSponge - @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = - dtγ * ( - ᶠtridiagonal_matrix_c3 ⋅ - DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - - DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) - ) - (I_u₃,) - else - @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = - dtγ * ᶠtridiagonal_matrix_c3 ⋅ - DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - (I_u₃,) - end - elseif p.atmos.rayleigh_sponge isa RayleighSponge - ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = - matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] - @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = - dtγ * -DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) - - (I_u₃,) - end - end - + A.cache.b_krylov .= b + ldiv!(A.cache.x_krylov, A, A.cache.b_krylov) + x .= A.cache.x_krylov end diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 20f437a911..5562303ee2 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -9,8 +9,8 @@ import ClimaAtmos as CA import LinearAlgebra import ClimaCore.Fields import ClimaTimeSteppers as CTS - import ClimaDiagnostics +import DiffEqBase: KeywordArgSilent function get_atmos(config::AtmosConfig, params) (; turbconv_params) = params @@ -369,56 +369,29 @@ function get_surface_setup(parsed_args) return getproperty(SurfaceConditions, Symbol(parsed_args["surface_setup"]))() end -is_explicit_CTS_algo_type(alg_or_tableau) = - alg_or_tableau <: CTS.ERKAlgorithmName - -is_imex_CTS_algo_type(alg_or_tableau) = - alg_or_tableau <: CTS.IMEXARKAlgorithmName - -is_implicit_type(alg_or_tableau) = is_imex_CTS_algo_type(alg_or_tableau) - -is_imex_CTS_algo(::CTS.IMEXAlgorithm) = true -is_imex_CTS_algo(::CTS.RosenbrockAlgorithm) = true -is_imex_CTS_algo(::SciMLBase.AbstractODEAlgorithm) = false - -is_implicit(ode_algo) = is_imex_CTS_algo(ode_algo) - -is_rosenbrock(::SciMLBase.AbstractODEAlgorithm) = false -use_transform(ode_algo) = - !(is_imex_CTS_algo(ode_algo) || is_rosenbrock(ode_algo)) - -additional_integrator_kwargs(::SciMLBase.AbstractODEAlgorithm) = (; - adaptive = false, - progress = isinteractive(), - progress_steps = isinteractive() ? 1 : 1000, -) -import DiffEqBase -additional_integrator_kwargs(::CTS.DistributedODEAlgorithm) = (; - kwargshandle = DiffEqBase.KeywordArgSilent, # allow custom kwargs - adjustfinal = true, - # TODO: enable progress bars in ClimaTimeSteppers -) - -is_cts_algo(::SciMLBase.AbstractODEAlgorithm) = false -is_cts_algo(::CTS.DistributedODEAlgorithm) = true - -function jac_kwargs(ode_algo, Y, atmos, parsed_args) - if is_implicit(ode_algo) - A = ImplicitEquationJacobian( - Y, - atmos; - approximate_solve_iters = parsed_args["approximate_linear_solve_iters"], - transform_flag = use_transform(ode_algo), - ) - if use_transform(ode_algo) - return (; jac_prototype = A, Wfact_t = Wfact!) +jac_kwargs(ode_algo, Y, p, parsed_args) = + if ode_algo isa Union{CTS.IMEXAlgorithm, CTS.RosenbrockAlgorithm} + use_exact_jacobian = parsed_args["use_exact_jacobian"] + always_update_exact_jacobian = + parsed_args["n_steps_update_exact_jacobian"] == 0 + approximate_solve_iters = parsed_args["approximate_linear_solve_iters"] + jacobian_algorithm = if parsed_args["debug_approximate_jacobian"] + DebugJacobian( + ApproxJacobian(; approximate_solve_iters); + use_exact_jacobian, + always_update_exact_jacobian, + ) else - return (; jac_prototype = A, Wfact = Wfact!) + use_exact_jacobian ? + ExactJacobian(; always_update_exact_jacobian) : + ApproxJacobian(; approximate_solve_iters) end + @info "Jacobian algorithm: $(dump_string(jacobian_algorithm))" + A = ImplicitEquationJacobian(jacobian_algorithm, Y, p) + (; jac_prototype = A, Wfact = update_jacobian!) else - return NamedTuple() + (;) end -end #= ode_configuration(Y, parsed_args) @@ -427,17 +400,14 @@ Returns the ode algorithm =# function ode_configuration(::Type{FT}, parsed_args) where {FT} ode_name = parsed_args["ode_algo"] - alg_or_tableau = getproperty(CTS, Symbol(ode_name)) - @info "Using ODE config: `$alg_or_tableau`" - if ode_name == "SSPKnoth" - return CTS.RosenbrockAlgorithm(CTS.tableau(CTS.SSPKnoth())) - end - - if is_explicit_CTS_algo_type(alg_or_tableau) - return CTS.ExplicitAlgorithm(alg_or_tableau()) - elseif !is_implicit_type(alg_or_tableau) - return alg_or_tableau() - elseif is_imex_CTS_algo_type(alg_or_tableau) + ode_algo_name = getproperty(CTS, Symbol(ode_name)) + @info "Using ODE config: `$ode_algo_name`" + return if ode_algo_name <: CTS.RosenbrockAlgorithmName + CTS.RosenbrockAlgorithm(CTS.tableau(ode_algo_name())) + elseif ode_algo_name <: CTS.ERKAlgorithmName + CTS.ExplicitAlgorithm(ode_algo_name()) + else + @assert ode_algo_name <: CTS.IMEXARKAlgorithmName newtons_method = CTS.NewtonsMethod(; max_iters = parsed_args["max_newton_iters_ode"], krylov_method = if parsed_args["use_krylov_method"] @@ -466,9 +436,7 @@ function ode_configuration(::Type{FT}, parsed_args) where {FT} nothing end, ) - return CTS.IMEXAlgorithm(alg_or_tableau(), newtons_method) - else - return alg_or_tableau(; linsolve = linsolve!) + CTS.IMEXAlgorithm(ode_algo_name(), newtons_method) end end @@ -522,29 +490,20 @@ function get_sim_info(config::AtmosConfig) end function args_integrator(parsed_args, Y, p, tspan, ode_algo, callback) - (; atmos, dt) = p + (; dt) = p dt_save_to_sol = time_to_seconds(parsed_args["dt_save_to_sol"]) s = @timed_str begin func = if parsed_args["split_ode"] - implicit_func = SciMLBase.ODEFunction( - implicit_tendency!; - jac_kwargs(ode_algo, Y, atmos, parsed_args)..., - tgrad = (∂Y∂t, Y, p, t) -> (∂Y∂t .= 0), - ) - if is_cts_algo(ode_algo) - CTS.ClimaODEFunction(; - T_exp_T_lim! = remaining_tendency!, - T_imp! = implicit_func, - # Can we just pass implicit_tendency! and jac_prototype etc.? - lim! = limiters_func!, - dss!, - post_explicit! = set_precomputed_quantities!, - post_implicit! = set_precomputed_quantities!, - ) - else - SciMLBase.SplitFunction(implicit_func, remaining_tendency!) - end + jacobian = jac_kwargs(ode_algo, Y, p, parsed_args) + CTS.ClimaODEFunction(; + T_exp_T_lim! = remaining_tendency!, + T_imp! = SciMLBase.ODEFunction(implicit_tendency!; jacobian...), + lim! = limiters_func!, + dss!, + post_explicit! = set_precomputed_quantities!, + post_implicit! = set_precomputed_quantities!, + ) # TODO: Separate the Wfact! and jac_prototype kwargs from T_imp!. else remaining_tendency! # should be total_tendency! end @@ -560,7 +519,13 @@ function args_integrator(parsed_args, Y, p, tspan, ode_algo, callback) end # ensure that tspan[2] is always saved @info "dt_save_to_sol: $dt_save_to_sol, length(saveat): $(length(saveat))" args = (problem, ode_algo) - kwargs = (; saveat, callback, dt, additional_integrator_kwargs(ode_algo)...) + kwargs = (; + saveat, + callback, + dt, + kwargshandle = KeywordArgSilent, # allow custom kwargs + adjustfinal = true, + ) return (args, kwargs) end @@ -719,6 +684,13 @@ function get_simulation(config::AtmosConfig) end @info "init integrator: $s" + if ( + config.parsed_args["use_exact_jacobian"] || + config.parsed_args["debug_approximate_jacobian"] + ) && config.parsed_args["n_steps_update_exact_jacobian"] == 0 + update_exact_jacobian!(integrator) + end # TODO: Is this necessary for initialization? + if config.parsed_args["enable_diagnostics"] s = @timed_str begin integrator = ClimaDiagnostics.IntegratorWithDiagnostics( diff --git a/src/solver/types.jl b/src/solver/types.jl index 81a6deacdf..1e54fa9062 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -175,6 +175,7 @@ struct MixingLength{FT} tke::FT buoy::FT end +MixingLength(args...) = MixingLength{promote_type(typeof.(args)...)}(args...) abstract type AbstractEDMF end diff --git a/src/utils/abbreviations.jl b/src/utils/abbreviations.jl index 8b686bc2e9..6b5dea700d 100644 --- a/src/utils/abbreviations.jl +++ b/src/utils/abbreviations.jl @@ -1,4 +1,5 @@ -using ClimaCore: Geometry, Operators, MatrixFields +using ClimaCore: Geometry, Operators +using ClimaCore.MatrixFields # Alternatively, we could use Vec₁₂₃, Vec³, etc., if that is more readable. const C1 = Geometry.Covariant1Vector diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl index c3ac7a4d4a..79f82fcb2b 100644 --- a/src/utils/utilities.jl +++ b/src/utils/utilities.jl @@ -4,6 +4,7 @@ import ClimaComms import ClimaCore: Spaces, Topologies, Fields, Geometry import LinearAlgebra: norm_sqr +import Statistics: quantile is_energy_var(symbol) = symbol in (:ρe_tot, :ρae_tot) is_momentum_var(symbol) = symbol in (:uₕ, :ρuₕ, :u₃, :ρw) @@ -283,28 +284,40 @@ end import Dates """ - time_and_units_str(x::Real) - -Returns a truncated string of time and units, -given a time `x` in Seconds. -""" -time_and_units_str(x::Real) = - trunc_time(string(compound_period(x, Dates.Second))) + time_and_units_str(seconds) + +Converts the given number of seconds into a string that contains a numerical +value (rounded to the nearest 2 decimal digits) and its corresponding units. +""" +function time_and_units_str(seconds) + # Convert seconds from Float32 to Float64 to avoid round-off errors. + precise_seconds = seconds isa Float32 ? Float64(seconds) : seconds + full_period = compound_period(precise_seconds) + isempty(Dates.periods(full_period)) && return "0 Seconds" + whole_period = Dates.periods(full_period)[1] + whole_period_value = Dates.value(whole_period) + remaining_period = full_period - whole_period + no_remaining_digits = isempty(Dates.periods(remaining_period)) + period_units = + string(typeof(whole_period).name.name) * + (whole_period_value == 1 && no_remaining_digits ? "" : "s") + no_remaining_digits && return "$whole_period_value $period_units" + remaining_period_value = + Dates.value(convert(Dates.Nanosecond, remaining_period)) / + Dates.value(convert(Dates.Nanosecond, typeof(whole_period)(1))) + remaining_period_value < 0.01 && return "$whole_period_value $period_units" + remaining_digits = lpad(round(Int, remaining_period_value * 100), 2, '0') + return "$whole_period_value.$remaining_digits $period_units" +end """ - compound_period(x::Real, ::Type{T}) where {T <: Dates.Period} + compound_period(seconds) -A canonicalized `Dates.CompoundPeriod` given a real value -`x`, and its units via the period type `T`. +A `Dates.CompoundPeriod` that represents the given number of seconds (rounded to +the nearest nanosecond). """ -function compound_period(x::Real, ::Type{T}) where {T <: Dates.Period} - nf = Dates.value(convert(Dates.Nanosecond, T(1))) - ns = Dates.Nanosecond(ceil(x * nf)) - return Dates.canonicalize(Dates.CompoundPeriod(ns)) -end - -trunc_time(s::String) = count(',', s) > 1 ? join(split(s, ",")[1:2], ",") : s - +compound_period(seconds) = + Dates.canonicalize(Dates.Nanosecond(round(Int, seconds * 10^9))) function prettymemory(b) if b < 1024 @@ -332,6 +345,19 @@ macro timed_str(ex) end end +""" + dump_string(x) + +Returns a string that contains the output of `dump(x)`. +""" +function dump_string(x) + buffer = IOBuffer() + dump(buffer, x) + result = String(take!(buffer)) + close(buffer) + return result +end + struct AllNothing end const all_nothing = AllNothing() Base.getproperty(::AllNothing, ::Symbol) = nothing @@ -352,6 +378,19 @@ function horizontal_integral_at_boundary(f::Fields.Field) sum(f ./ Fields.Δz_field(axes(f)) .* 2) # TODO: is there a way to ensure this is derived from face z? The 2d topology doesn't contain this info end +""" + smooth_maximum(array, [percentile_rank]) + +Filters out all non-negligible values in the given `array` and computes the +specified percentile of the result. +""" +function smooth_maximum(array, percentile_rank = 1 - 0.1^ndims(array)) + minimum_filtered_value = maximum(abs, array) / 1e6 + filtered_array_values = filter(>(minimum_filtered_value), map(abs, array)) + isempty(filtered_array_values) && return zero(eltype(array)) + return quantile(filtered_array_values, percentile_rank) +end + """ gaussian_smooth(arr, sigma = 1) diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl index 721bc78b0e..9a2dcd8036 100644 --- a/src/utils/variable_manipulations.jl +++ b/src/utils/variable_manipulations.jl @@ -62,9 +62,9 @@ function as (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2`. """ sgs_weight_function(a, a_half) = - if a < 0 + if a <= 0 # autodiff generates NaNs when a is 0 zero(a) - elseif a > 1 + elseif a > min(1, 42 * a_half) # autodiff generates NaNs when a is large one(a) else (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2 diff --git a/test/dependencies.jl b/test/dependencies.jl index 92e507ee52..267fac7c68 100644 --- a/test/dependencies.jl +++ b/test/dependencies.jl @@ -38,6 +38,7 @@ known_dependencies = Set([ "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", + "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", @@ -48,6 +49,7 @@ known_dependencies = Set([ "Dates", "DiffEqBase", "FastGaussQuadrature", + "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts",