move long running tilted v into its own test, run on cluster

.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:

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"

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

test/Land/Model/test_overland_flow.jl → 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