diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 03b0e0bf8a2..fc5d2ce924c 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: @@ -2012,6 +2021,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/docs/list_of_apis.jl b/docs/list_of_apis.jl index afd70d1a25f..d2fd7830cc8 100644 --- a/docs/list_of_apis.jl +++ b/docs/list_of_apis.jl @@ -22,6 +22,9 @@ apis = [ "APIs/Land/SoilWaterParameterizations.md", "Soil Heat Parameterizations" => "APIs/Land/SoilHeatParameterizations.md", + "Surface Flow" => "APIs/Land/SurfaceFlow.md", + "Radiative Boundary Conditions" => + "APIs/Land/RadiativeEnergyFlux.md", ], "Common" => [ "Orientations" => "APIs/Common/Orientations.md", diff --git a/docs/src/APIs/Land/SurfaceFlow.md b/docs/src/APIs/Land/SurfaceFlow.md new file mode 100644 index 00000000000..78e470663e3 --- /dev/null +++ b/docs/src/APIs/Land/SurfaceFlow.md @@ -0,0 +1,17 @@ +# SurfaceFlow + +```@meta +CurrentModule = ClimateMachine.Land.SurfaceFlow +``` + +## SurfaceFlow types and functions +```@docs +OverlandFlowModel +NoSurfaceFlowModel +surface_boundary_flux! +surface_boundary_state! +calculate_velocity +Precip +VolumeAdvection +SurfaceWaterArea +``` \ No newline at end of file diff --git a/src/Driver/driver_configs.jl b/src/Driver/driver_configs.jl index dc7e27d66e6..40c1961ae32 100644 --- a/src/Driver/driver_configs.jl +++ b/src/Driver/driver_configs.jl @@ -670,7 +670,8 @@ function MultiColumnLandModel( zmin = zero(FT), array_type = ClimateMachine.array_type(), mpicomm = MPI.COMM_WORLD, - boundary = ((3, 3), (3, 3), (1, 2)), + boundary = ((3, 4), (5, 6), (1, 2)), + solver_type = ExplicitSolverType(), periodicity = (false, false, false), meshwarp = (x...) -> identity(x), numerical_flux_first_order = CentralNumericalFluxFirstOrder(), diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 2e8283d119b..ab25e771cee 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -30,13 +30,14 @@ import ..BalanceLaws: compute_gradient_flux!, nodal_init_state_auxiliary!, init_state_prognostic!, - nodal_update_auxiliary_state! + nodal_update_auxiliary_state!, + wavespeed using ..DGMethods: LocalGeometry, DGModel export LandModel """ - LandModel{PS, S, LBC, SRC, IS} <: BalanceLaw + LandModel{PS, S, SF, LBC, SRC, SRCDT, IS} <: BalanceLaw A BalanceLaw for land modeling. Users may over-ride prescribed default values for each field. @@ -46,19 +47,23 @@ Users may over-ride prescribed default values for each field. LandModel( param_set, soil; + surface, boundary_conditions, source, + source_dt, init_state_prognostic ) # Fields $(DocStringExtensions.FIELDS) """ -struct LandModel{PS, S, LBC, SRC, SRCDT, IS} <: BalanceLaw +struct LandModel{PS, S, SF, LBC, SRC, SRCDT, IS} <: BalanceLaw "Parameter set" param_set::PS "Soil model" soil::S + "Surface Flow model" + surface::SF "struct of boundary conditions" boundary_conditions::LBC "Source Terms (Problem specific source terms)" @@ -75,6 +80,7 @@ parameter_set(m::LandModel) = m.param_set LandModel( param_set::AbstractParameterSet, soil::BalanceLaw; + surface::BalanceLaw = NoSurfaceFlowModel(), boundary_conditions::LBC = (), source::SRC = (), init_state_prognostic::IS = nothing @@ -85,6 +91,7 @@ Constructor for the LandModel structure. function LandModel( param_set::AbstractParameterSet, soil::BalanceLaw; + surface::BalanceLaw = NoSurfaceFlowModel(), boundary_conditions::LBC = LandDomainBC(), source::SRC = (), init_state_prognostic::IS = nothing, @@ -94,6 +101,7 @@ function LandModel( land = ( param_set, soil, + surface, boundary_conditions, source, source_dt, @@ -106,26 +114,32 @@ end function vars_state(land::LandModel, st::Prognostic, FT) @vars begin soil::vars_state(land.soil, st, FT) + surface::vars_state(land.surface, st, FT) end end function vars_state(land::LandModel, st::Auxiliary, FT) @vars begin + x::FT + y::FT z::FT soil::vars_state(land.soil, st, FT) + surface::vars_state(land.surface, st, FT) end end function vars_state(land::LandModel, st::Gradient, FT) @vars begin soil::vars_state(land.soil, st, FT) + surface::vars_state(land.surface, st, FT) end end function vars_state(land::LandModel, st::GradientFlux, FT) @vars begin soil::vars_state(land.soil, st, FT) + surface::vars_state(land.surface, st, FT) end end @@ -135,8 +149,11 @@ function nodal_init_state_auxiliary!( tmp::Vars, geom::LocalGeometry, ) + aux.x = geom.coord[1] + aux.y = geom.coord[2] aux.z = geom.coord[3] land_init_aux!(land, land.soil, aux, geom) + land_init_aux!(land, land.surface, aux, geom) end function flux_first_order!( @@ -145,8 +162,19 @@ function flux_first_order!( state::Vars, aux::Vars, t::Real, - directions, -) end + direction, +) + tend = Flux{FirstOrder}() + _args = (; state, aux, t, direction) + args = merge(_args, (precomputed = precompute(land, _args, tend),)) + + map(prognostic_vars(land)) do prog + var, name = get_prog_state(flux, prog) + val = Σfluxes(prog, eq_tends(prog, land, tend), land, args) + setproperty!(var, name, val) + end + nothing +end function compute_gradient_argument!( @@ -210,6 +238,7 @@ function nodal_update_auxiliary_state!( t::Real, ) land_nodal_update_auxiliary_state!(land, land.soil, state, aux, t) + land_nodal_update_auxiliary_state!(land, land.surface, state, aux, t) end @@ -247,7 +276,6 @@ function init_state_prognostic!( end include("prog_types.jl") - include("RadiativeEnergyFlux.jl") using .RadiativeEnergyFlux include("SoilWaterParameterizations.jl") @@ -260,11 +288,19 @@ include("soil_heat.jl") include("Runoff.jl") using .Runoff include("land_bc.jl") +include("SurfaceFlow.jl") +using .SurfaceFlow include("soil_bc.jl") - include("prognostic_vars.jl") - -include("source.jl") include("land_tendencies.jl") +include("source.jl") + +function wavespeed(m::LandModel, n⁻, state::Vars, aux::Vars, t::Real, direction) + FT = eltype(state) + h = max(state.surface.height, FT(0.0)) + v = calculate_velocity(m.surface, aux.x, aux.y, h) + speed = norm(v) + return speed +end end # Module diff --git a/src/Land/Model/SurfaceFlow.jl b/src/Land/Model/SurfaceFlow.jl new file mode 100644 index 00000000000..d6a21d50aa2 --- /dev/null +++ b/src/Land/Model/SurfaceFlow.jl @@ -0,0 +1,247 @@ +module SurfaceFlow + +using DocStringExtensions +using UnPack +using ..Land +using ..VariableTemplates +using ..BalanceLaws +import ..BalanceLaws: + BalanceLaw, + prognostic_vars, + flux, + source, + precompute, + eq_tends, + vars_state, + Prognostic, + Auxiliary, + Gradient, + GradientFlux + +using ...DGMethods: LocalGeometry, DGModel +using StaticArrays: SVector + +export OverlandFlowModel, + NoSurfaceFlowModel, + surface_boundary_flux!, + surface_boundary_state!, + calculate_velocity, + Precip, + VolumeAdvection, + SurfaceWaterHeight + +""" + SurfaceWaterHeight <: AbstractPrognosticVariable + +The prognostic variable type for the 2d overland flow model. Used only for +dispatching on. +""" +struct SurfaceWaterHeight <: AbstractPrognosticVariable end + + +""" + NoSurfaceFlowModel <: BalanceLaw + +The default surface flow model, which does not add any prognostic variables +to the land model and therefore does not model surface flow. +""" +struct NoSurfaceFlowModel <: BalanceLaw end + +""" + OverlandFlowModel{Sx,Sy,M} <: BalanceLaw + +The 2D overland flow model, with a prognostic variable +equal to the height of the surface water. + +This model simulates the depth-averaged shallow water equation under the +kinematic approximation, and employs Manning's relationship to relate +velocity to the height of the water. +# Fields +$(DocStringExtensions.FIELDS) +""" +struct OverlandFlowModel{Sx, Sy, M} <: BalanceLaw + "Slope in x direction; field(x,y), unitless" + slope_x::Sx + "Slope in y direction; field(x,y), unitless" + slope_y::Sy + "Mannings coefficient; field(x,y), units of s/m^(1/3)" + mannings::M +end + +function OverlandFlowModel( + slope_x::Function, + slope_y::Function; + mannings::Function = (x, y) -> convert(eltype(x), 0.03), +) + args = (slope_x, slope_y, mannings) + return OverlandFlowModel{typeof.(args)...}(args...) +end + + +""" + calculate_velocity(surface, x::Real, y::Real, h::Real) + +Given the surface flow model, calculate the velocity of the flow +based on height, slope, and Manning's coefficient. +""" +function calculate_velocity(surface, x::Real, y::Real, h::Real) + FT = eltype(h) + sx = FT(surface.slope_x(x, y)) + sy = FT(surface.slope_y(x, y)) + mannings_coeff = FT(surface.mannings(x, y)) + coeff = h^FT(2 / 3) / mannings_coeff + #The velocity direction is opposite the slope vector (∂_x z, ∂_y z) + return SVector( + -sign(sx) * coeff * sqrt(abs(sx)), + -sign(sy) * coeff * sqrt(abs(sy)), + zero(FT), + ) +end + +vars_state(surface::OverlandFlowModel, st::Prognostic, FT) = @vars(height::FT) + +function Land.land_init_aux!( + land::LandModel, + surface::Union{NoSurfaceFlowModel, OverlandFlowModel}, + aux, + geom::LocalGeometry, +) end + +function Land.land_nodal_update_auxiliary_state!( + land::LandModel, + surface::Union{NoSurfaceFlowModel, OverlandFlowModel}, + state, + aux, + t, +) end + +""" + VolumeAdvection <: TendencyDef{Flux{FirstOrder}} + +A first order flux type for the overland flow model. +""" +struct VolumeAdvection <: TendencyDef{Flux{FirstOrder}} end + + +""" + flux(::SurfaceWaterHeight, ::VolumeAdvection, land::LandModel, args,) + +A first order flux method for the OverlandFlow model, adding in advection of water volume. +""" +function flux(::SurfaceWaterHeight, ::VolumeAdvection, land::LandModel, args) + @unpack state, aux = args + x = aux.x + y = aux.y + height = max(state.surface.height, eltype(aux)(0.0)) + v = calculate_velocity(land.surface, x, y, height) + return height * v +end + +# Boundary Conditions + +# General case - to be used with bc::NoBC +function surface_boundary_flux!( + nf, + bc::Land.AbstractBoundaryConditions, + m::Union{NoSurfaceFlowModel, OverlandFlowModel}, + land::LandModel, + _..., +) end + +function surface_boundary_state!( + nf, + bc::Land.AbstractBoundaryConditions, + m::Union{NoSurfaceFlowModel, OverlandFlowModel}, + land::LandModel, + _..., +) end + +""" + surface_boundary_flux!( + nf, + bc::Land.Dirichlet, + model::OverlandFlowModel, + land::LandModel, + _..., + ) + +The surface boundary flux function for the OverlandFlow model, which +does nothing if Dirichlet conditions are chosen in order to not +overconstrain the solution. +""" +function surface_boundary_flux!( + nf, + bc::Land.Dirichlet, + model::OverlandFlowModel, + land::LandModel, + _..., +) end + +""" + surface_boundary_state!( + nf, + bc::Land.Dirichlet, + model::OverlandFlowModel, + land::LandModel, + state⁺::Vars, + aux⁺::Vars, + nM, + state⁻::Vars, + aux⁻::Vars, + t, + _..., + ) + +The surface boundary state function for the OverlandFlow model when +Dirichlet conditions are chosen. This should be equivalent to +outflow boundary conditions, also referred to as gradient outlet +conditions. +""" +function surface_boundary_state!( + nf, + bc::Land.Dirichlet, + model::OverlandFlowModel, + land::LandModel, + state⁺::Vars, + aux⁺::Vars, + nM, + state⁻::Vars, + aux⁻::Vars, + t, + _..., +) + bc_function = bc.state_bc + state⁺.surface.height = bc_function(aux⁻, t) +end + +""" + Precip <: TendencyDef{Source} + +A source term for overland flow where a prescribed net precipitation +drives the overland flow. +""" +struct Precip{FT, F} <: TendencyDef{Source} + precip::F + + function Precip{FT}(precip::F) where {FT, F} + new{FT, F}(precip) + end +end + +function (p::Precip{FT})(x, y, t) where {FT} + FT(p.precip(x, y, t)) +end + + +prognostic_vars(::Precip) = (SurfaceWaterHeight(),) + + +precompute(source_type::Precip, land::LandModel, args, tt::Source) = + NamedTuple() + +function source(::SurfaceWaterHeight, s::Precip, land::LandModel, args) + @unpack aux, t = args + return s(aux.x, aux.y, t) +end + +end diff --git a/src/Land/Model/land_bc.jl b/src/Land/Model/land_bc.jl index 5190ec4062e..530f5b53e91 100644 --- a/src/Land/Model/land_bc.jl +++ b/src/Land/Model/land_bc.jl @@ -141,17 +141,20 @@ boundary conditions. The user should supply an instance of `LandComponentBC` for each piece, as needed. If none is supplied, the default for is `NoBC` for each subcomponent. At a minimum, both top and bottom boundary conditions should be supplied. Whether or not to include the lateral faces depends on the configuration of the domain. """ -Base.@kwdef struct LandDomainBC{TBC, BBC, LBC} +Base.@kwdef struct LandDomainBC{TBC, BBC, MinXBC, MaxXBC, MinYBC, MaxYBC} "surface boundary conditions" surface_bc::TBC = LandComponentBC() "bottom boundary conditions" bottom_bc::BBC = LandComponentBC() "lateral boundary conditions" - lateral_bc::LBC = LandComponentBC() + minx_bc::MinXBC = LandComponentBC() + maxx_bc::MaxXBC = LandComponentBC() + miny_bc::MinYBC = LandComponentBC() + maxy_bc::MaxYBC = LandComponentBC() end """ - LandComponentBC{SW, SH} + LandComponentBC{SW, SH, SF} An object that holds the boundary conditions for each of the subcomponents of the land model. @@ -163,9 +166,11 @@ needs to define the BC for the components they wish to model. Base.@kwdef struct LandComponentBC{ SW <: AbstractBoundaryConditions, SH <: AbstractBoundaryConditions, + SF <: AbstractBoundaryConditions, } soil_water::SW = NoBC() soil_heat::SH = NoBC() + surface::SF = NoBC() end """ @@ -177,8 +182,15 @@ faces, as defined in the Driver configuration. """ function boundary_conditions(land::LandModel) bc = land.boundary_conditions - mytuple = (bc.bottom_bc, bc.surface_bc, bc.lateral_bc) - # faces labeled integer 1,2,3 are bottom, top, lateral sides. + mytuple = ( + bc.bottom_bc, + bc.surface_bc, + bc.minx_bc, + bc.maxx_bc, + bc.miny_bc, + bc.maxy_bc, + ) + # faces labeled integer 1,2 are bottom, top, lateral sides are 3, 4, 5, 6 return mytuple end @@ -211,6 +223,7 @@ end function land_boundary_state!(nf, bc::LandComponentBC, land, args...) soil_boundary_state!(nf, bc.soil_water, land.soil.water, land, args...) soil_boundary_state!(nf, bc.soil_heat, land.soil.heat, land, args...) + surface_boundary_state!(nf, bc.surface, land.surface, land, args...) end function boundary_state!( @@ -248,4 +261,5 @@ end function land_boundary_flux!(nf, bc::LandComponentBC, land, args...) soil_boundary_flux!(nf, bc.soil_water, land.soil.water, land, args...) soil_boundary_flux!(nf, bc.soil_heat, land.soil.heat, land, args...) + surface_boundary_flux!(nf, bc.surface, land.surface, land, args...) end diff --git a/src/Land/Model/land_tendencies.jl b/src/Land/Model/land_tendencies.jl index cd4b0892dc8..059ffdb5a5b 100644 --- a/src/Land/Model/land_tendencies.jl +++ b/src/Land/Model/land_tendencies.jl @@ -9,19 +9,34 @@ eq_tends(pv::AbstractPrognosticVariable, m::LandModel, tt::Source) = ##### First order fluxes ##### -eq_tends(::PV, ::LandModel, ::Flux{FirstOrder}) where {PV} = () +eq_tends(pv::PV, land::LandModel, tt::Flux{FirstOrder}) where {PV} = ( + eq_tends(pv, land.soil.heat, tt)..., + eq_tends(pv, land.soil.water, tt)..., + eq_tends(pv, land.surface, tt)..., +) +eq_tends(::PV, ::AbstractSoilComponentModel, ::Flux{FirstOrder}) where {PV} = () +eq_tends(::PV, ::NoSurfaceFlowModel, ::Flux{FirstOrder}) where {PV} = () + +eq_tends(::SurfaceWaterHeight, ::OverlandFlowModel, ::Flux{FirstOrder}) = + (VolumeAdvection(),) ##### ##### Second order fluxes ##### # Empty by default + eq_tends(pv::PV, ::AbstractSoilComponentModel, ::Flux{SecondOrder}) where {PV} = () +eq_tends(pv::PV, ::OverlandFlowModel, ::Flux{SecondOrder}) where {PV} = () +eq_tends(pv::PV, ::NoSurfaceFlowModel, ::Flux{SecondOrder}) where {PV} = () -eq_tends(pv::PV, land::LandModel, tt::Flux{SecondOrder}) where {PV} = - (eq_tends(pv, land.soil.heat, tt)..., eq_tends(pv, land.soil.water, tt)...) +eq_tends(pv::PV, land::LandModel, tt::Flux{SecondOrder}) where {PV} = ( + eq_tends(pv, land.soil.heat, tt)..., + eq_tends(pv, land.soil.water, tt)..., + eq_tends(pv, land.surface, tt)..., +) eq_tends( pv::PV, @@ -30,7 +45,6 @@ eq_tends( ) where {PV <: VolumetricInternalEnergy} = (eq_tends(pv, land.soil.heat, land.soil.water, tt)...,) -# TODO: move to soi_heat.jl? eq_tends( ::VolumetricInternalEnergy, ::SoilHeatModel, @@ -38,7 +52,6 @@ eq_tends( ::Flux{SecondOrder}, ) = (DiffHeatFlux(), DarcyDrivenHeatFlux()) -# TODO: move to soi_water.jl? eq_tends( ::VolumetricInternalEnergy, ::SoilHeatModel, diff --git a/src/Land/Model/prognostic_vars.jl b/src/Land/Model/prognostic_vars.jl index acd841e4cf8..4cbd9e8f98d 100644 --- a/src/Land/Model/prognostic_vars.jl +++ b/src/Land/Model/prognostic_vars.jl @@ -7,10 +7,18 @@ prognostic_vars(water::SoilWaterModel) = (VolumetricLiquidFraction(), VolumetricIceFraction()) prognostic_vars(heat::PrescribedTemperatureModel) = () prognostic_vars(heat::SoilHeatModel) = (VolumetricInternalEnergy(),) +prognostic_vars(surface::NoSurfaceFlowModel) = () +prognostic_vars(surface::OverlandFlowModel) = (SurfaceWaterHeight(),) + +prognostic_vars(land::LandModel) = ( + prognostic_vars(land.soil.water)..., + prognostic_vars(land.soil.heat)..., + prognostic_vars(land.surface)..., +) + -prognostic_vars(land::LandModel) = - (prognostic_vars(land.soil.water)..., prognostic_vars(land.soil.heat)...) get_prog_state(state, ::VolumetricLiquidFraction) = (state.soil.water, :ϑ_l) get_prog_state(state, ::VolumetricIceFraction) = (state.soil.water, :θ_i) get_prog_state(state, ::VolumetricInternalEnergy) = (state.soil.heat, :ρe_int) +get_prog_state(state, ::SurfaceWaterHeight) = (state.surface, :height) diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index b4e83e9bcbe..9968f799e9d 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -32,7 +32,7 @@ end function precompute(source_type::PhaseChange, land::LandModel, args, tt::Source) @unpack state, diffusive, aux, t, direction = args - + FT = eltype(state) param_set = parameter_set(land) diff --git a/test/Land/Model/Artifacts.toml b/test/Land/Model/Artifacts.toml index e449680d943..d3f6b6c8eeb 100644 --- a/test/Land/Model/Artifacts.toml +++ b/test/Land/Model/Artifacts.toml @@ -3,3 +3,6 @@ git-tree-sha1 = "ff73fa6a0b6a807e71a6921f7ef7d0befe776edd" [richards_sand] git-tree-sha1 = "b0dc82dd02159c646e909bfb61170d3b9dc347f3" + +[tiltedv] +git-tree-sha1 = "db27235cb7ce2b7674607876da15d1635906b512" \ No newline at end of file diff --git a/test/Land/Model/prescribed_twice.jl b/test/Land/Model/prescribed_twice.jl index fd0ff045ff2..3d38a6d3162 100644 --- a/test/Land/Model/prescribed_twice.jl +++ b/test/Land/Model/prescribed_twice.jl @@ -91,8 +91,8 @@ using ClimateMachine.BalanceLaws: aux, vars_state(m, Auxiliary(), FT), ) - #Make sure it runs, and that there are no state variables, and only "z" as aux. + #Make sure it runs, and that there are no state variables, and only "x,y,z" as aux. @test t == timeend - @test size(Base.collect(keys(aux_vars)))[1] == 1 + @test size(Base.collect(keys(aux_vars)))[1] == 3 @test size(Base.collect(keys(state_vars)))[1] == 0 end diff --git a/test/Land/Model/runtests.jl b/test/Land/Model/runtests.jl index 46754a96f1b..c7e204e5ffb 100644 --- a/test/Land/Model/runtests.jl +++ b/test/Land/Model/runtests.jl @@ -4,5 +4,6 @@ using Test, Pkg include("test_water_parameterizations.jl") include("prescribed_twice.jl") include("freeze_thaw_alone.jl") + include("test_overland_flow_analytic.jl") include("test_physical_bc.jl") end diff --git a/test/Land/Model/test_bc_3d.jl b/test/Land/Model/test_bc_3d.jl index 159da38c37b..0464af9a85d 100644 --- a/test/Land/Model/test_bc_3d.jl +++ b/test/Land/Model/test_bc_3d.jl @@ -44,14 +44,16 @@ using ClimateMachine.BalanceLaws: bottom_flux = (aux, t) -> bottom_flux_amplitude * sin(f * t) * aux.soil.water.K surface_state = (aux, t) -> eltype(aux)(0.2) + lateral_state = (aux, t) -> eltype(aux)(0.0) ϑ_l0 = (aux) -> eltype(aux)(0.2) bc = LandDomainBC( bottom_bc = LandComponentBC(soil_water = Neumann(bottom_flux)), surface_bc = LandComponentBC(soil_water = Dirichlet(surface_state)), - lateral_bc = LandComponentBC( - soil_water = Neumann((aux, t) -> eltype(aux)(0.0)), - ), + minx_bc = LandComponentBC(soil_water = Neumann(lateral_state)), + maxx_bc = LandComponentBC(soil_water = Neumann(lateral_state)), + miny_bc = LandComponentBC(soil_water = Neumann(lateral_state)), + maxy_bc = LandComponentBC(soil_water = Neumann(lateral_state)), ) soil_water_model = SoilWaterModel(FT; initialϑ_l = ϑ_l0) soil_heat_model = PrescribedTemperatureModel() diff --git a/test/Land/Model/test_overland_flow_analytic.jl b/test/Land/Model/test_overland_flow_analytic.jl new file mode 100644 index 00000000000..0c1943eb05c --- /dev/null +++ b/test/Land/Model/test_overland_flow_analytic.jl @@ -0,0 +1,288 @@ +using MPI +using OrderedCollections +using StaticArrays +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.MPIStateArrays +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 + + +# 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() + 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) + m_surface = NoSurfaceFlowModel() + + sources = () + m = LandModel( + param_set, + m_soil; + surface = m_surface, + source = sources, + init_state_prognostic = init_land_model!, + ) + + N_poly = 5 + nelem_vert = 10 + + # Specify the domain boundaries + zmax = FT(0) + zmin = FT(-1) + + driver_config = ClimateMachine.SingleStackConfiguration( + "LandModel", + N_poly, + nelem_vert, + zmax, + param_set, + m; + zmin = zmin, + numerical_flux_first_order = CentralNumericalFluxFirstOrder(), + ) + + t0 = FT(0) + timeend = FT(10) + dt = FT(1) + + solver_config = ClimateMachine.SolverConfiguration( + t0, + timeend, + driver_config, + ode_dt = dt, + ) + mygrid = solver_config.dg.grid + Q = solver_config.Q + aux = solver_config.dg.state_auxiliary + + ClimateMachine.invoke!(solver_config;) + + t = ODESolvers.gettime(solver_config.solver) + state_vars = SingleStackUtils.get_vars_from_nodal_stack( + mygrid, + Q, + vars_state(m, Prognostic(), FT), + ) + aux_vars = SingleStackUtils.get_vars_from_nodal_stack( + mygrid, + aux, + vars_state(m, Auxiliary(), FT), + ) + #Make sure it runs, and that there are no state variables, and only "x,y,z" as aux. + @test t == timeend + @test size(Base.collect(keys(aux_vars)))[1] == 3 + @test size(Base.collect(keys(state_vars)))[1] == 0 +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 + function warp_constant_slope( + xin, + yin, + zin; + topo_max = 0.2, + zmin = -0.1, + xmax = 400, + ) + FT = eltype(xin) + zmax = FT((FT(1.0) - xin / xmax) * topo_max) + alpha = FT(1.0) - zmax / zmin + zout = zmin + (zin - zmin) * alpha + 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) + + m_surface = OverlandFlowModel( + (x, y) -> eltype(x)(-0.0016), + (x, y) -> eltype(x)(0.0); + mannings = (x, y) -> eltype(x)(0.025), + ) + + + bc = LandDomainBC( + minx_bc = LandComponentBC( + surface = Dirichlet((aux, t) -> eltype(aux)(0)), + ), + ) + + function init_land_model!(land, state, aux, localgeo, time) + state.surface.height = eltype(state)(0) + end + + # units in m / s + precip(x, y, t) = t < (30 * 60) ? 1.4e-5 : 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 = 1 + xres = FT(2.286) + yres = FT(0.25) + zres = FT(0.1) + # Specify the domain boundaries. + zmax = FT(0) + zmin = FT(-0.1) + xmax = FT(182.88) + ymax = FT(1.0) + topo_max = FT(0.0016 * xmax) + + driver_config = ClimateMachine.MultiColumnLandModel( + "LandModel", + (N_poly, N_poly), + (xres, yres, zres), + xmax, + ymax, + zmax, + param_set, + m; + zmin = zmin, + meshwarp = (x...) -> warp_constant_slope( + x...; + topo_max = topo_max, + zmin = zmin, + xmax = xmax, + ), + ) + + t0 = FT(0) + timeend = FT(60 * 60) + dt = FT(10) + + solver_config = ClimateMachine.SolverConfiguration( + t0, + timeend, + driver_config, + ode_dt = dt, + ) + mygrid = solver_config.dg.grid + Q = solver_config.Q + + h_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :height) + 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) + h = Q[:, h_index, :] + all_vars = Dict{String, Array}("t" => [t], "h" => h) + dons[iostep[1]] = all_vars + iostep[1] += 1 + return + end + + ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) + + # Compare flowrate analytical derivation + aux = solver_config.dg.state_auxiliary + + # get all nodal points at the max X bound of the domain + mask = Array(aux[:, 1, :] .== 182.88) + n_outputs = length(dons) + # get prognostic variable height from nodal state (m^2) + height = [mean(Array(dons[k]["h"])[mask[:]]) for k in 1:n_outputs] + # get similation timesteps (s) + time_data = [dons[l]["t"][1] for l in 1:n_outputs] + + + alpha = sqrt(0.0016) / 0.025 + i = 1.4e-5 + L = xmax + m = 5 / 3 + t_c = (L * i^(1 - m) / alpha)^(1 / m) + t_r = 30 * 60 + + q = height .^ (m) .* alpha + + function g(m, y, i, t_r, L, alpha, t) + output = L / alpha - y^(m) / i - y^(m - 1) * m * (t - t_r) + return output + end + + function dg(m, y, i, t_r, L, alpha, t) + output = -y^(m - 1) * m / i - y^(m - 2) * m * (m - 1) * (t - t_r) + return output + end + + function analytic(t, alpha, t_c, t_r, i, L, m) + if t < t_c + return alpha * (i * t)^(m) + end + + if t <= t_r && t > t_c + return alpha * (i * t_c)^(m) + end + + if t > t_r + yL = (i * (t - t_r)) + delta = 1 + error = g(m, yL, i, t_r, L, alpha, t) + while abs(error) > 1e-4 + delta = + -g(m, yL, i, t_r, L, alpha, t) / + dg(m, yL, i, t_r, L, alpha, t) + yL = yL + delta + error = g(m, yL, i, t_r, L, alpha, t) + end + return alpha * yL^m + + end + + end + + 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 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..575a83e8746 --- /dev/null +++ b/test/Land/Model/test_overland_flow_vcatchment.jl @@ -0,0 +1,186 @@ +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) + + 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; 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.height = eltype(state)(0) + end + + # units in m / s + precip(x, y, t) = t < eltype(t)(90 * 60) ? eltype(t)(3e-6) : eltype(t)(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(), + ) + + 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 + + h_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :height) + 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) + h = Q[:, h_index, :] + all_vars = Dict{String, Array}("t" => [t], "h" => h) + dons[iostep[1]] = all_vars + iostep[1] += 1 + return + end + + ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) + + aux = solver_config.dg.state_auxiliary + x = Array(aux[:, 1, :]) + y = Array(aux[:, 2, :]) + z = Array(aux[:, 3, :]) + # Get points at outlet (y = ymax) + mask2 = (Float64.(y .== 1000.0)) .== 1 + n_outputs = length(dons) + function compute_Q(h, xv) + height = max.(h, 0.0) + 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]["h"])[:][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 < 5e-7 +end