diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 8ba637d08cf..78518be8920 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -762,6 +762,15 @@ steps: queue: central slurm_ntasks: 1 + - label: "cpu_land_overland_flow_vcatchment" + key: "cpu_land_overland_flow_vcatchment" + command: + - "mpiexec julia --color=yes --project test/Land/Model/test_overland_flow_vcatchment.jl" + agents: + config: cpu + queue: central + slurm_ntasks: 1 + - label: "gpu_thermodynamics" key: "gpu_thermodynamics" command: @@ -1992,6 +2001,16 @@ steps: slurm_ntasks: 1 slurm_gres: "gpu:1" + - label: "gpu_land_overland_flow_vcatchment" + key: "gpu_land_overland_flow_vcatchment" + command: + - "mpiexec julia --color=yes --project test/Land/Model/test_overland_flow_vcatchment.jl" + agents: + config: gpu + queue: central + slurm_ntasks: 1 + slurm_gres: "gpu:1" + - label: "gpu_unittests" key: "gpu_unittests" command: diff --git a/test/Land/Model/Artifacts.toml b/test/Land/Model/Artifacts.toml index 8a685a2256a..94aebbb28dc 100644 --- a/test/Land/Model/Artifacts.toml +++ b/test/Land/Model/Artifacts.toml @@ -1,9 +1,8 @@ [richards] git-tree-sha1 = "ff73fa6a0b6a807e71a6921f7ef7d0befe776edd" -[tiltedv] -git-tree-sha1 = "db27235cb7ce2b7674607876da15d1635906b512" - [richards_sand] git-tree-sha1 = "b0dc82dd02159c646e909bfb61170d3b9dc347f3" +[tiltedv] +git-tree-sha1 = "db27235cb7ce2b7674607876da15d1635906b512" diff --git a/test/Land/Model/runtests.jl b/test/Land/Model/runtests.jl index 803ba0f3ed2..c7e204e5ffb 100644 --- a/test/Land/Model/runtests.jl +++ b/test/Land/Model/runtests.jl @@ -4,6 +4,6 @@ using Test, Pkg include("test_water_parameterizations.jl") include("prescribed_twice.jl") include("freeze_thaw_alone.jl") - include("test_overland_flow.jl") + include("test_overland_flow_analytic.jl") include("test_physical_bc.jl") end diff --git a/test/Land/Model/test_overland_flow.jl b/test/Land/Model/test_overland_flow_analytic.jl similarity index 57% rename from test/Land/Model/test_overland_flow.jl rename to test/Land/Model/test_overland_flow_analytic.jl index fb76164a8ae..aea9ee39477 100644 --- a/test/Land/Model/test_overland_flow.jl +++ b/test/Land/Model/test_overland_flow_analytic.jl @@ -27,6 +27,8 @@ using ClimateMachine.BalanceLaws: BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state using ClimateMachine.ArtifactWrappers + +# Test that the land model with no surface flow works correctly @testset "NoSurfaceFlow Model" begin function init_land_model!(land, state, aux, localgeo, time) end ClimateMachine.init() @@ -67,7 +69,7 @@ using ClimateMachine.ArtifactWrappers ) t0 = FT(0) - timeend = FT(60) + timeend = FT(10) dt = FT(1) solver_config = ClimateMachine.SolverConfiguration( @@ -100,12 +102,9 @@ using ClimateMachine.ArtifactWrappers end +# Constant slope analytical test case defined as Model 1 / Eqn 6 +# DOI: 10.1061/(ASCE)0733-9429(2007)133:2(217) @testset "Analytical Overland Model" begin - # Analytical test case defined as Model 1 in DOI: - # 10.1061/(ASCE)0733-9429(2007)133:2(217) - # Eqn 6 - - function warp_constant_slope( xin, yin, @@ -283,195 +282,9 @@ end end - # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a GPU array) + q = Array(q) # copy to host if GPU array @test sqrt_rmse_over_max_q = sqrt(mean( (analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q) .^ 2.0, )) / maximum(q) < 3e-3 end - - -@testset "V Catchment Maxwell River Model" begin - tv_dataset = ArtifactWrapper( - @__DIR__, - isempty(get(ENV, "CI", "")), - "tiltedv", - ArtifactFile[ArtifactFile( - url = "https://caltech.box.com/shared/static/qi2gftjw2vu2j66b0tyfef427xxj3ug7.csv", - filename = "TiltedVOutput.csv", - ),], - ) - tv_dataset_path = get_data_folder(tv_dataset) - - function warp_tilted_v(xin, yin, zin) - FT = eltype(xin) - slope_sides = FT(0.05) - slope_v = (0.02) - zbase = slope_v * (yin - FT(1000)) - zleft = FT(0.0) - zright = FT(0.0) - if xin < FT(800) - zleft = slope_sides * (xin - FT(800)) - end - if xin > FT(820) - zright = slope_sides * (FT(1) + (xin - FT(820)) / FT(800)) - end - zout = zbase + zleft + zright - x, y, z = xin, yin, zout - return x, y, z - end - ClimateMachine.init() - FT = Float64 - - soil_water_model = PrescribedWaterModel() - soil_heat_model = PrescribedTemperatureModel() - soil_param_functions = nothing - - m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) - - function x_slope(x, y) - MFT = eltype(x) - if x < MFT(800) - MFT(-0.05) - elseif x <= MFT(820) - MFT(0) - else - MFT(0.05) - end - end - - function y_slope(x, y) - MFT = eltype(x) - MFT(-0.02) - end - - function channel_mannings(x, y) - MFT = eltype(x) - return x >= MFT(800) && x <= MFT(820) ? MFT(2.5 * 60 * 10^-3) : - MFT(2.5 * 60 * 10^-4) - end - - m_surface = OverlandFlowModel( - x_slope, - y_slope, - (x, y) -> eltype(x)(1); - mannings = channel_mannings, - ) - - bc = LandDomainBC( - miny_bc = LandComponentBC( - surface = Dirichlet((aux, t) -> eltype(aux)(0)), - ), - minx_bc = LandComponentBC( - surface = Dirichlet((aux, t) -> eltype(aux)(0)), - ), - maxx_bc = LandComponentBC( - surface = Dirichlet((aux, t) -> eltype(aux)(0)), - ), - ) - - function init_land_model!(land, state, aux, localgeo, time) - state.surface.area = eltype(state)(0) - end - - # units in m / s - precip(x, y, t) = t < (90 * 60) ? 3e-6 : 0.0 - - sources = (Precip{FT}(precip),) - - m = LandModel( - param_set, - m_soil; - surface = m_surface, - boundary_conditions = bc, - source = sources, - init_state_prognostic = init_land_model!, - ) - - N_poly_hori = 1 - N_poly_vert = 1 - xres = FT(20) - yres = FT(20) - zres = FT(1) - # Specify the domain boundaries. - zmax = FT(1) - zmin = FT(0) - xmax = FT(1620) - ymax = FT(1000) - - - driver_config = ClimateMachine.MultiColumnLandModel( - "LandModel", - (N_poly_hori, N_poly_vert), - (xres, yres, zres), - xmax, - ymax, - zmax, - param_set, - m; - zmin = zmin, - numerical_flux_first_order = RusanovNumericalFlux(), - # meshwarp = (x...) -> warp_tilted_v(x...), - ) - - t0 = FT(0) - timeend = FT(180 * 60) - dt = FT(0.5) - - solver_config = ClimateMachine.SolverConfiguration( - t0, - timeend, - driver_config, - ode_dt = dt, - ) - mygrid = solver_config.dg.grid - Q = solver_config.Q - - area_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :area) - n_outputs = 60 - - every_x_simulation_time = ceil(Int, timeend / n_outputs) - - dons = Dict([k => Dict() for k in 1:n_outputs]...) - - iostep = [1] - callback = GenericCallbacks.EveryXSimulationTime( - every_x_simulation_time, - ) do (init = false) - t = ODESolvers.gettime(solver_config.solver) - area = Q[:, area_index, :] - all_vars = Dict{String, Array}("t" => [t], "area" => area) - dons[iostep[1]] = all_vars - iostep[1] += 1 - return - end - - ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) - - aux = solver_config.dg.state_auxiliary - x = aux[:, 1, :] - y = aux[:, 2, :] - z = aux[:, 3, :] - # Get points at outlet (y = ymax) - mask2 = (Float64.(aux[:, 2, :] .== 1000.0)) .== 1 - n_outputs = length(dons) - function compute_Q(a, xv) - height = max.(a, 0.0) ./ 1.0 # the width = 1m here - i.e. we are solving actually for h, not area - v = calculate_velocity(m_surface, xv, 1000.0, height)# put in y = 1000.0 - speed = sqrt(v[1]^2.0 + v[2]^2.0 + v[3]^2.0) - Q_outlet = speed .* height .* 60.0 # multiply by 60 so it is per minute, not per second - return Q_outlet - end - # We divide by 4 because we have 4 nodal points with the same value at each x (z = 0, 1) - # Multiply by xres because the solution at each point roughly represents that the solution for that range in x - Q = - [ - sum(compute_Q.(Array(dons[k]["area"])[:][mask2[:]], x[mask2[:]])) - for k in 1:n_outputs - ] ./ 4.0 .* xres - data = joinpath(tv_dataset_path, "TiltedVOutput.csv") - ds_tv = readdlm(data, ',') - error = sqrt(mean(Q .- ds_tv)) - @test error < 1e-10 - -end diff --git a/test/Land/Model/test_overland_flow_vcatchment.jl b/test/Land/Model/test_overland_flow_vcatchment.jl new file mode 100644 index 00000000000..c6b8a841571 --- /dev/null +++ b/test/Land/Model/test_overland_flow_vcatchment.jl @@ -0,0 +1,211 @@ +using Test +using Statistics +using DelimitedFiles + +using CLIMAParameters +struct EarthParameterSet <: AbstractEarthParameterSet end +const param_set = EarthParameterSet() + +using ClimateMachine +using ClimateMachine.Land +using ClimateMachine.Land.SurfaceFlow +using ClimateMachine.Land.SoilWaterParameterizations +using ClimateMachine.Mesh.Topologies +using ClimateMachine.Mesh.Grids +using ClimateMachine.DGMethods +using ClimateMachine.DGMethods.NumericalFluxes +using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry +using ClimateMachine.GenericCallbacks +using ClimateMachine.ODESolvers +using ClimateMachine.VariableTemplates +using ClimateMachine.SingleStackUtils +using ClimateMachine.BalanceLaws: + BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state +using ClimateMachine.ArtifactWrappers + + +@testset "V Catchment Maxwell River Model" begin + tv_dataset = ArtifactWrapper( + @__DIR__, + isempty(get(ENV, "CI", "")), + "tiltedv", + ArtifactFile[ArtifactFile( + url = "https://caltech.box.com/shared/static/qi2gftjw2vu2j66b0tyfef427xxj3ug7.csv", + filename = "TiltedVOutput.csv", + ),], + ) + tv_dataset_path = get_data_folder(tv_dataset) + + function warp_tilted_v(xin, yin, zin) + FT = eltype(xin) + slope_sides = FT(0.05) + slope_v = (0.02) + zbase = slope_v * (yin - FT(1000)) + zleft = FT(0.0) + zright = FT(0.0) + if xin < FT(800) + zleft = slope_sides * (xin - FT(800)) + end + if xin > FT(820) + zright = slope_sides * (FT(1) + (xin - FT(820)) / FT(800)) + end + zout = zbase + zleft + zright + x, y, z = xin, yin, zout + return x, y, z + end + + ClimateMachine.init() + FT = Float64 + + soil_water_model = PrescribedWaterModel() + soil_heat_model = PrescribedTemperatureModel() + soil_param_functions = nothing + + m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) + + function x_slope(x, y) + MFT = eltype(x) + if x < MFT(800) + MFT(-0.05) + elseif x <= MFT(820) + MFT(0) + else + MFT(0.05) + end + end + + function y_slope(x, y) + MFT = eltype(x) + MFT(-0.02) + end + + function channel_mannings(x, y) + MFT = eltype(x) + return x >= MFT(800) && x <= MFT(820) ? MFT(2.5 * 60 * 10^-3) : + MFT(2.5 * 60 * 10^-4) + end + + m_surface = OverlandFlowModel( + x_slope, + y_slope, + (x, y) -> eltype(x)(1); + mannings = channel_mannings, + ) + + bc = LandDomainBC( + miny_bc = LandComponentBC( + surface = Dirichlet((aux, t) -> eltype(aux)(0)), + ), + minx_bc = LandComponentBC( + surface = Dirichlet((aux, t) -> eltype(aux)(0)), + ), + maxx_bc = LandComponentBC( + surface = Dirichlet((aux, t) -> eltype(aux)(0)), + ), + ) + + function init_land_model!(land, state, aux, localgeo, time) + state.surface.area = eltype(state)(0) + end + + # units in m / s + precip(x, y, t) = t < (90 * 60) ? 3e-6 : 0.0 + + sources = (Precip{FT}(precip),) + + m = LandModel( + param_set, + m_soil; + surface = m_surface, + boundary_conditions = bc, + source = sources, + init_state_prognostic = init_land_model!, + ) + + N_poly_hori = 1 + N_poly_vert = 1 + xres = FT(20) + yres = FT(20) + zres = FT(1) + # Specify the domain boundaries. + zmax = FT(1) + zmin = FT(0) + xmax = FT(1620) + ymax = FT(1000) + + + driver_config = ClimateMachine.MultiColumnLandModel( + "LandModel", + (N_poly_hori, N_poly_vert), + (xres, yres, zres), + xmax, + ymax, + zmax, + param_set, + m; + zmin = zmin, + numerical_flux_first_order = RusanovNumericalFlux(), + # meshwarp = (x...) -> warp_tilted_v(x...), + ) + + t0 = FT(0) + timeend = FT(180 * 60) + dt = FT(0.5) + + solver_config = ClimateMachine.SolverConfiguration( + t0, + timeend, + driver_config, + ode_dt = dt, + ) + mygrid = solver_config.dg.grid + Q = solver_config.Q + + area_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :area) + n_outputs = 60 + + every_x_simulation_time = ceil(Int, timeend / n_outputs) + + dons = Dict([k => Dict() for k in 1:n_outputs]...) + + iostep = [1] + callback = GenericCallbacks.EveryXSimulationTime( + every_x_simulation_time, + ) do (init = false) + t = ODESolvers.gettime(solver_config.solver) + area = Q[:, area_index, :] + all_vars = Dict{String, Array}("t" => [t], "area" => area) + dons[iostep[1]] = all_vars + iostep[1] += 1 + return + end + + ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) + + aux = solver_config.dg.state_auxiliary + x = aux[:, 1, :] + y = aux[:, 2, :] + z = aux[:, 3, :] + # Get points at outlet (y = ymax) + mask2 = (Float64.(aux[:, 2, :] .== 1000.0)) .== 1 + n_outputs = length(dons) + function compute_Q(a, xv) + height = max.(a, 0.0) ./ 1.0 # the width = 1m here - i.e. we are solving actually for h, not area + v = calculate_velocity(m_surface, xv, 1000.0, height)# put in y = 1000.0 + speed = sqrt(v[1]^2.0 + v[2]^2.0 + v[3]^2.0) + Q_outlet = speed .* height .* 60.0 # multiply by 60 so it is per minute, not per second + return Q_outlet + end + # We divide by 4 because we have 4 nodal points with the same value at each x (z = 0, 1) + # Multiply by xres because the solution at each point roughly represents that the solution for that range in x + Q = + [ + sum(compute_Q.(Array(dons[k]["area"])[:][mask2[:]], x[mask2[:]])) + for k in 1:n_outputs + ] ./ 4.0 .* xres + data = joinpath(tv_dataset_path, "TiltedVOutput.csv") + ds_tv = readdlm(data, ',') + error = sqrt(mean(Q .- ds_tv)) + @test error < 1e-10 + +end