diff --git a/docs/list_diagnostics.jl b/docs/list_diagnostics.jl index 61ce19fe1a..4d76cb4f6f 100644 --- a/docs/list_diagnostics.jl +++ b/docs/list_diagnostics.jl @@ -2,4 +2,4 @@ diagnostics = [ "Available diagnostics" => "diagnostics/available_diagnostics.md", "For developers" => "diagnostics/developers_diagnostics.md", "For users" => "diagnostics/users_diagnostics.md", -] +] diff --git a/experiments/standalone/Bucket/global_bucket_diagnostics.jl b/experiments/standalone/Bucket/global_bucket_diagnostics.jl deleted file mode 100644 index 9f0c06d9ae..0000000000 --- a/experiments/standalone/Bucket/global_bucket_diagnostics.jl +++ /dev/null @@ -1,212 +0,0 @@ -import SciMLBase -import ClimaComms -@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends -using CairoMakie -using Dates -using DelimitedFiles -using Statistics - -import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput - -import ClimaTimeSteppers as CTS -import NCDatasets -using ClimaCore -using ClimaCore: Remapping, Geometry -import ClimaParams as CP -import ClimaComms -import ClimaLand -import ClimaLand.Parameters as LP -using ClimaLand.Bucket: - BucketModel, BucketModelParameters, PrescribedSurfaceAlbedo -using ClimaLand.Domains: coordinates, Column -using ClimaLand: - initialize, - make_update_aux, - make_exp_tendency, - make_set_initial_cache, - PrescribedAtmosphere, - PrescribedRadiativeFluxes - -PROFILING = false -try - import Profile, ProfileCanvas - global PROFILING = true - @info "ProfileCanvas found, running with profiler" -catch -end - -""" - compute_extrema(v) - -Computes and returns the minimum value in `v` and -the maximum value in `v`, as a tuple, assuming that -`v` is a vector of arrays. -""" -function compute_extrema(v) - maxes = [maximum(u) for u in v] - mins = [minimum(u) for u in v] - return (minimum(mins), maximum(maxes)) -end - -anim_plots = false -FT = Float64; -context = ClimaComms.context() -earth_param_set = LP.LandParameters(FT); -outdir = joinpath( - pkgdir(ClimaLand), - "experiments/standalone/Bucket/artifacts_temporalmap", -) -device_suffix = - typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? - "cpu" : "gpu" -!ispath(outdir) && mkpath(outdir) -# Use separate output directory for CPU and GPU runs to avoid race condition -device_suffix = - typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? - "cpu" : "gpu" -t0 = 0.0; -tf = 2 * 86400; -Δt = 3600.0; - - -# function setup_prob(t0, tf, Δt) -# We set up the problem in a function so that we can make multiple copies (for profiling) - -# Set up simulation domain -soil_depth = FT(3.5) -bucket_domain = ClimaLand.Domains.SphericalShell(; - radius = FT(6.3781e6), - depth = soil_depth, - nelements = (10, 10), # this failed with (50,10) - npolynomial = 1, - dz_tuple = FT.((1.0, 0.05)), -) -ref_time = DateTime(2005) - -# Initialize parameters -σS_c = FT(0.2) -W_f = FT(0.15) -z_0m = FT(1e-2) -z_0b = FT(1e-3) -κ_soil = FT(0.7) -ρc_soil = FT(2e6) -τc = FT(3600) - -surface_space = bucket_domain.space.surface -# Construct albedo parameter object using temporal map -albedo = PrescribedSurfaceAlbedo{FT}(ref_time, t0, surface_space) - -bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc) - -# Precipitation: -precip = (t) -> 0 -snow_precip = (t) -> -5e-7 * (t < 1 * 86400) -# Diurnal temperature variations: -T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2) -# Constant otherwise: -u_atmos = (t) -> 3.0 -q_atmos = (t) -> 0.001 -h_atmos = FT(2) -P_atmos = (t) -> 101325 - -bucket_atmos = PrescribedAtmosphere( - TimeVaryingInput(precip), - TimeVaryingInput(snow_precip), - TimeVaryingInput(T_atmos), - TimeVaryingInput(u_atmos), - TimeVaryingInput(q_atmos), - TimeVaryingInput(P_atmos), - ref_time, - h_atmos, - earth_param_set, -) - -# Prescribed radiation -- a prescribed downwelling SW diurnal cycle, with a -# peak at local noon, and a prescribed downwelling LW radiative -# flux, assuming the air temperature is on average 275 degrees -# K with a diurnal amplitude of 5 degrees K: -SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0) -LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4 -bucket_rad = PrescribedRadiativeFluxes( - FT, - TimeVaryingInput(SW_d), - TimeVaryingInput(LW_d), - ref_time, -) - - -model = BucketModel( - parameters = bucket_parameters, - domain = bucket_domain, - atmosphere = bucket_atmos, - radiation = bucket_rad, -) - -Y, p, _coords = initialize(model) - -Y.bucket.T .= FT(270) -Y.bucket.W .= FT(0.05) -Y.bucket.Ws .= FT(0.0) -Y.bucket.σS .= FT(0.08) - -set_initial_cache! = make_set_initial_cache(model) -set_initial_cache!(p, Y, t0) -exp_tendency! = make_exp_tendency(model) -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction((T_exp!) = exp_tendency!, (dss!) = ClimaLand.dss!), - Y, - (t0, tf), - p, -) -saveat = collect(t0:(Δt * 3):tf) -saved_values = (; - t = Array{Float64}(undef, length(saveat)), - saveval = Array{NamedTuple}(undef, length(saveat)), -) -saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat) -updateat = copy(saveat) -updatefunc = ClimaLand.make_update_drivers(bucket_atmos, bucket_rad) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb) - -# return prob, cb, saveat, saved_values -# end - -# prob, cb, saveat, saved_values = setup_prob(t0, tf, Δt); -timestepper = CTS.RK4() -ode_algo = CTS.ExplicitAlgorithm(timestepper) - -#### ClimaDiagnostics #### -using ClimaDiagnostics - -if isdir("output") - rm("output", force = true, recursive = true) -end - -mkdir("output") -output_dir = "output/" - -space = bucket_domain.space.subsurface - -nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir) - -diags = ClimaLand.CLD.default_diagnostics(model, 1.0; output_writer = nc_writer) - -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) - -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) - -sol_test = SciMLBase.solve(prob, ode_algo; dt = Δt, callback = diag_cb) - -#### ClimaAnalysis #### - -using ClimaAnalysis -simdir = ClimaAnalysis.SimDir("output") -println(summary(simdir)) -alpha = get(simdir; short_name = "alpha") -alpha.dims -import ClimaAnalysis.Visualize as viz -fig = CairoMakie.Figure(size = (400, 600)) -viz.plot!(fig, alpha, time = 2 * 86400) -CairoMakie.save("alpha.png", fig) diff --git a/test/diagnostics/diagnostics_tests.jl b/test/diagnostics/diagnostics_tests.jl new file mode 100644 index 0000000000..1486e2c1f9 --- /dev/null +++ b/test/diagnostics/diagnostics_tests.jl @@ -0,0 +1,212 @@ +using Test +using ClimaDiagnostics +using ClimaAnalysis +using ClimaUtilities +import SciMLBase +import ClimaComms +@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends +using CairoMakie +using Dates +using DelimitedFiles +using Statistics + +import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput + +import ClimaTimeSteppers as CTS +import NCDatasets +using ClimaCore +using ClimaCore: Remapping, Geometry +import ClimaParams as CP +import ClimaComms +import ClimaLand +import ClimaLand.Parameters as LP +using ClimaLand.Bucket: + BucketModel, BucketModelParameters, PrescribedSurfaceAlbedo +using ClimaLand.Domains: coordinates, Column +using ClimaLand: + initialize, + make_update_aux, + make_exp_tendency, + make_set_initial_cache, + PrescribedAtmosphere, + PrescribedRadiativeFluxes + +@testset "Bucket model global run" begin + + PROFILING = false + try + import Profile, ProfileCanvas + global PROFILING = true + @info "ProfileCanvas found, running with profiler" + catch + end + + """ + compute_extrema(v) + + Computes and returns the minimum value in `v` and + the maximum value in `v`, as a tuple, assuming that + `v` is a vector of arrays. + """ + function compute_extrema(v) + maxes = [maximum(u) for u in v] + mins = [minimum(u) for u in v] + return (minimum(mins), maximum(maxes)) + end + + anim_plots = false + FT = Float64; + context = ClimaComms.context() + earth_param_set = LP.LandParameters(FT); + outdir = joinpath( + pkgdir(ClimaLand), + "experiments/standalone/Bucket/artifacts_temporalmap", + ) + device_suffix = + typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? + "cpu" : "gpu" + !ispath(outdir) && mkpath(outdir) + # Use separate output directory for CPU and GPU runs to avoid race condition + device_suffix = + typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? + "cpu" : "gpu" + t0 = 0.0; + tf = 2 * 86400; + Δt = 3600.0; + + + # function setup_prob(t0, tf, Δt) + # We set up the problem in a function so that we can make multiple copies (for profiling) + + # Set up simulation domain + soil_depth = FT(3.5) + bucket_domain = ClimaLand.Domains.SphericalShell(; + radius = FT(6.3781e6), + depth = soil_depth, + nelements = (10, 10), # this failed with (50,10) + npolynomial = 1, + dz_tuple = FT.((1.0, 0.05)), + ) + ref_time = DateTime(2005) + + # Initialize parameters + σS_c = FT(0.2) + W_f = FT(0.15) + z_0m = FT(1e-2) + z_0b = FT(1e-3) + κ_soil = FT(0.7) + ρc_soil = FT(2e6) + τc = FT(3600) + + surface_space = bucket_domain.space.surface + # Construct albedo parameter object using temporal map + albedo = PrescribedSurfaceAlbedo{FT}(ref_time, t0, surface_space) + + bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc) + + # Precipitation: + precip = (t) -> 0 + snow_precip = (t) -> -5e-7 * (t < 1 * 86400) + # Diurnal temperature variations: + T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2) + # Constant otherwise: + u_atmos = (t) -> 3.0 + q_atmos = (t) -> 0.001 + h_atmos = FT(2) + P_atmos = (t) -> 101325 + + bucket_atmos = PrescribedAtmosphere( + TimeVaryingInput(precip), + TimeVaryingInput(snow_precip), + TimeVaryingInput(T_atmos), + TimeVaryingInput(u_atmos), + TimeVaryingInput(q_atmos), + TimeVaryingInput(P_atmos), + ref_time, + h_atmos, + earth_param_set, + ) + + # Prescribed radiation -- a prescribed downwelling SW diurnal cycle, with a + # peak at local noon, and a prescribed downwelling LW radiative + # flux, assuming the air temperature is on average 275 degrees + # K with a diurnal amplitude of 5 degrees K: + SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0) + LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4 + bucket_rad = PrescribedRadiativeFluxes( + FT, + TimeVaryingInput(SW_d), + TimeVaryingInput(LW_d), + ref_time, + ) + + + model = BucketModel( + parameters = bucket_parameters, + domain = bucket_domain, + atmosphere = bucket_atmos, + radiation = bucket_rad, + ) + + Y, p, _coords = initialize(model) + + Y.bucket.T .= FT(270) + Y.bucket.W .= FT(0.05) + Y.bucket.Ws .= FT(0.0) + Y.bucket.σS .= FT(0.08) + + set_initial_cache! = make_set_initial_cache(model) + set_initial_cache!(p, Y, t0) + exp_tendency! = make_exp_tendency(model) + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction((T_exp!) = exp_tendency!, (dss!) = ClimaLand.dss!), + Y, + (t0, tf), + p, + ) + saveat = collect(t0:(Δt * 3):tf) + saved_values = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), + ) + saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat) + updateat = copy(saveat) + updatefunc = ClimaLand.make_update_drivers(bucket_atmos, bucket_rad) + driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) + cb = SciMLBase.CallbackSet(driver_cb, saving_cb) + + # return prob, cb, saveat, saved_values + # end + + # prob, cb, saveat, saved_values = setup_prob(t0, tf, Δt); + timestepper = CTS.RK4() + ode_algo = CTS.ExplicitAlgorithm(timestepper) + + #### ClimaDiagnostics #### + generate_output_path("output"; style = ActiveLinkStyle, context = context) + output_dir = "output/" + + space = bucket_domain.space.subsurface + + nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir) + + diags = ClimaLand.CLD.default_diagnostics(model, 1.0; output_writer = nc_writer) + + diagnostic_handler = + ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) + + diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) + + sol_test = SciMLBase.solve(prob, ode_algo; dt = Δt, callback = diag_cb) + + #### ClimaAnalysis #### + + simdir = ClimaAnalysis.SimDir("output") + println(summary(simdir)) + alpha = get(simdir; short_name = "alpha") + alpha.dims + import ClimaAnalysis.Visualize as viz + fig = CairoMakie.Figure(size = (400, 600)) + viz.plot!(fig, alpha, time = 2 * 86400) + CairoMakie.save("alpha.png", fig) +end