From 9367e8a05fbd94df3fa5a41492ec25c12a2f2c3f Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Mon, 11 Jan 2021 11:19:34 -0700 Subject: [PATCH 01/24] initial land model with no-op river model --- src/Land/Model/LandModel.jl | 18 +++++- src/Land/Model/River.jl | 45 ++++++++++++++ test/Land/Model/runtests.jl | 1 + test/Land/Model/test_no_river.jl | 102 +++++++++++++++++++++++++++++++ 4 files changed, 163 insertions(+), 3 deletions(-) create mode 100644 src/Land/Model/River.jl create mode 100644 test/Land/Model/test_no_river.jl diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 048db296ab9..f2e5ed51c22 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -46,11 +46,13 @@ Users may over-ride prescribed default values for each field. # Fields $(DocStringExtensions.FIELDS) """ -struct LandModel{PS, S, LBC, SRC, IS} <: BalanceLaw +struct LandModel{PS, S, R, LBC, SRC, IS} <: BalanceLaw "Parameter set" param_set::PS "Soil model" soil::S + "River model" + river::R "struct of boundary conditions" boundary_conditions::LBC "Source Terms (Problem specific source terms)" @@ -72,13 +74,14 @@ Constructor for the LandModel structure. """ function LandModel( param_set::AbstractParameterSet, - soil::BalanceLaw; + soil::BalanceLaw, + river::BalanceLaw; boundary_conditions::LBC = LandDomainBC(), source::SRC = (), init_state_prognostic::IS = nothing, ) where {SRC, IS, LBC} @assert init_state_prognostic ≠ nothing - land = (param_set, soil, boundary_conditions, source, init_state_prognostic) + land = (param_set, soil, river, boundary_conditions, source, init_state_prognostic) return LandModel{typeof.(land)...}(land...) end @@ -86,6 +89,7 @@ end function vars_state(land::LandModel, st::Prognostic, FT) @vars begin soil::vars_state(land.soil, st, FT) + river::vars_state(land.river, st, FT) end end @@ -94,18 +98,21 @@ function vars_state(land::LandModel, st::Auxiliary, FT) @vars begin z::FT soil::vars_state(land.soil, st, FT) + river::vars_state(land.river, st, FT) end end function vars_state(land::LandModel, st::Gradient, FT) @vars begin soil::vars_state(land.soil, st, FT) + river::vars_state(land.river, st, FT) end end function vars_state(land::LandModel, st::GradientFlux, FT) @vars begin soil::vars_state(land.soil, st, FT) + river::vars_state(land.river, st, FT) end end @@ -117,6 +124,7 @@ function nodal_init_state_auxiliary!( ) aux.z = geom.coord[3] land_init_aux!(land, land.soil, aux, geom) + land_init_aux!(land, land.river, aux, geom) end function flux_first_order!( @@ -138,6 +146,7 @@ function compute_gradient_argument!( ) compute_gradient_argument!(land, land.soil, transform, state, aux, t) + compute_gradient_argument!(land, land.river, transform, state, aux, t) end function compute_gradient_flux!( @@ -190,6 +199,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.river, state, aux, t) end @@ -229,4 +239,6 @@ include("soil_water.jl") include("soil_heat.jl") include("soil_bc.jl") include("source.jl") +include("River.jl") +using .River end # Module diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl new file mode 100644 index 00000000000..b598bebafda --- /dev/null +++ b/src/Land/Model/River.jl @@ -0,0 +1,45 @@ +module River + +using ..VariableTemplates +using DocStringExtensions +using ..Land +using ..BalanceLaws +import ..BalanceLaws: + BalanceLaw, + vars_state, + flux_first_order!, + flux_second_order!, + source!, + boundary_conditions, + boundary_state!, + compute_gradient_argument!, + compute_gradient_flux!, + nodal_init_state_auxiliary!, + init_state_prognostic!, + nodal_update_auxiliary_state! + +export RiverModel, NoRiverModel +#abstract type AbstractRiverModel end + +struct NoRiverModel <: BalanceLaw end +struct RiverModel <: BalanceLaw end + +vars_state(water::RiverModel, st::Prognostic, FT) = @vars() +vars_state(water::RiverModel, st::Auxiliary, FT) = @vars() +vars_state(water::RiverModel, st::Gradient, FT) = @vars() +vars_state(water::RiverModel, st::GradientFlux, FT) = @vars() + + +function Land.land_init_aux!(land, river::BalanceLaw, aux, geom) +end + +function Land.compute_gradient_argument!(land, river::BalanceLaw, transform, state, aux, t) +end + +function Land.land_nodal_update_auxiliary_state!(land, river::BalanceLaw, state, aux, t) +end + + + + +end \ No newline at end of file diff --git a/test/Land/Model/runtests.jl b/test/Land/Model/runtests.jl index 1cb7b54f1b1..e532b5321eb 100644 --- a/test/Land/Model/runtests.jl +++ b/test/Land/Model/runtests.jl @@ -5,4 +5,5 @@ using Test, Pkg include("prescribed_twice.jl") include("freeze_thaw_alone.jl") include("test_runoff_functions.jl") + include("test_no_river.jl") end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl new file mode 100644 index 00000000000..3792d36c67e --- /dev/null +++ b/test/Land/Model/test_no_river.jl @@ -0,0 +1,102 @@ +# Test that the land model still runs, even with the lowest/simplest +# version of soil (prescribed heat and prescribed water - no state +# variables) +using MPI +using OrderedCollections +using StaticArrays +using Test + +using CLIMAParameters +struct EarthParameterSet <: AbstractEarthParameterSet end +const param_set = EarthParameterSet() + +using ClimateMachine +using ClimateMachine.Land +using ClimateMachine.Land.River +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 + +@testset "NoRiver Model" begin + ClimateMachine.init() + FT = Float64 + + function init_land_model!(land, state, aux, localgeo, time) end + + 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_river = NoRiverModel() + + sources = () + m = LandModel( + param_set, + m_soil, + m_river; + 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(60) + 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 "z" as aux. + @test t == timeend + @test size(Base.collect(keys(aux_vars)))[1] == 1 + @test size(Base.collect(keys(state_vars)))[1] == 0 +end From af669631e5e775795f54ab4e97a359d06fb73cd8 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Mon, 11 Jan 2021 12:06:58 -0700 Subject: [PATCH 02/24] wip model1 --- src/Land/Model/LandModel.jl | 4 +++- src/Land/Model/River.jl | 38 +++++++++++++++++++++++++++++++------ 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index f2e5ed51c22..3dbd724a5bc 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -134,7 +134,9 @@ function flux_first_order!( aux::Vars, t::Real, directions, -) end +) + flux_first_order!(land, land.river, flux, state, aux, t, directions) +end function compute_gradient_argument!( diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index b598bebafda..6e979ccef4e 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -22,24 +22,50 @@ export RiverModel, NoRiverModel #abstract type AbstractRiverModel end struct NoRiverModel <: BalanceLaw end -struct RiverModel <: BalanceLaw end -vars_state(water::RiverModel, st::Prognostic, FT) = @vars() +struct RiverModel{M} <: BalanceLaw + mannings::M +end + +function RiverModel(mannings) +end + +mannings_coeff(river::River Model) + +function calculate_velocity(river, x::Real, y::Real, h::Real) + sx = river.slope_x(x, y) + sy = river.slope_y(x, y) + magnitude = h^(2/3) / (river.mannings(x, y)) * sqrt(river.mag_slope(x, y)) + return SVector(sx*magnitude, sy*magnitude, 0) +end + +vars_state(water::RiverModel, st::Prognostic, FT) = @vars(height::FT) vars_state(water::RiverModel, st::Auxiliary, FT) = @vars() vars_state(water::RiverModel, st::Gradient, FT) = @vars() vars_state(water::RiverModel, st::GradientFlux, FT) = @vars() - -function Land.land_init_aux!(land, river::BalanceLaw, aux, geom) +function Land.land_init_aux!(land::LandModel, river::BalanceLaw, aux, geom) end -function Land.compute_gradient_argument!(land, river::BalanceLaw, transform, state, aux, t) +function Land.compute_gradient_argument!(land::LandModel, river::BalanceLaw, transform::Grad, state, aux, t) + # v = calculate_velocity(river, state.river.height,river) + # transform.river.grad_Q = state.river.height*v#*river.width(aux.x,aux.y) end -function Land.land_nodal_update_auxiliary_state!(land, river::BalanceLaw, state, aux, t) + +function Land.land_nodal_update_auxiliary_state!(land::LandModel, river::BalanceLaw, state, aux, t) end +function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) + x = aux.x + y = aux.y + width = river.width(x,y) + height = state.river.area / width + v = calculate_velocity(river, x, y ,height) + Q = state.river.area * v + flux.river.height = Q +end end \ No newline at end of file From 7f13a1a2ab8ccc3d6e506a517a6a2cfe2393b716 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Mon, 11 Jan 2021 12:30:09 -0700 Subject: [PATCH 03/24] wip model1 --- .gitignore | 1 + src/Land/Model/River.jl | 38 ++++++++++++++++++++++++++------------ 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/.gitignore b/.gitignore index e3b3fb30abe..1a5e0f0f722 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,7 @@ *.jl.*.cov *.jl.mem *.DS_Store +*.vscode *~ TAGS diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 6e979ccef4e..6571f1f58a5 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -19,24 +19,39 @@ import ..BalanceLaws: nodal_update_auxiliary_state! export RiverModel, NoRiverModel -#abstract type AbstractRiverModel end struct NoRiverModel <: BalanceLaw end -struct RiverModel{M} <: BalanceLaw +struct RiverModel{M,Sx,Sy,MS,W} <: BalanceLaw mannings::M + slope_x::Sx + slope_y::Sy + mag_slope::MS + width::W end -function RiverModel(mannings) +function RiverModel( + slope_x::Function, + slope_y::Function, + mag_slope::Function, + width::Function; + mannings::Function = (x, y) -> convert(eltype(x), 0.03)) + args = ( + slope_x, + slope_y, + mag_slope, + width, + mannings + ) + return RiverModel{typeof.(args)...}(args...) end -mannings_coeff(river::River Model) - function calculate_velocity(river, x::Real, y::Real, h::Real) - sx = river.slope_x(x, y) - sy = river.slope_y(x, y) - magnitude = h^(2/3) / (river.mannings(x, y)) * sqrt(river.mag_slope(x, y)) - return SVector(sx*magnitude, sy*magnitude, 0) + FT = eltype(h) + sx = FT(river.slope_x(x, y)) + sy = FT(river.slope_y(x, y)) + magnitude = h^FT(2/3) / (river.mannings(x, y) * sqrt(river.mag_slope(x, y))) + return SVector(sx * magnitude, sy * magnitude, zero(FT)) end vars_state(water::RiverModel, st::Prognostic, FT) = @vars(height::FT) @@ -58,9 +73,8 @@ end function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) - x = aux.x - y = aux.y - width = river.width(x,y) + x, y = aux.x, aux.y + width = river.width(x, y) height = state.river.area / width v = calculate_velocity(river, x, y ,height) Q = state.river.area * v From 40cafc7cfe2da48155d42949bd5a756de5804a2d Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Wed, 13 Jan 2021 16:04:53 -0700 Subject: [PATCH 04/24] wip bc --- src/Driver/driver_configs.jl | 2 +- src/Land/Model/River.jl | 2 +- src/Land/Model/land_bc.jl | 11 +++- test/Land/Model/test_no_river.jl | 103 +++++++++++++++++++++++++++++++ 4 files changed, 113 insertions(+), 5 deletions(-) diff --git a/src/Driver/driver_configs.jl b/src/Driver/driver_configs.jl index ecdc9bd5dfd..6fcbabb83c0 100644 --- a/src/Driver/driver_configs.jl +++ b/src/Driver/driver_configs.jl @@ -491,7 +491,7 @@ 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), diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 6571f1f58a5..0d83bd78e96 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -30,7 +30,7 @@ struct RiverModel{M,Sx,Sy,MS,W} <: BalanceLaw width::W end -function RiverModel( +function axyerModel( slope_x::Function, slope_y::Function, mag_slope::Function, diff --git a/src/Land/Model/land_bc.jl b/src/Land/Model/land_bc.jl index 42455dedcd4..e3f68f08f82 100644 --- a/src/Land/Model/land_bc.jl +++ b/src/Land/Model/land_bc.jl @@ -112,7 +112,10 @@ Base.@kwdef struct LandDomainBC{TBC, BBC, LBC} "bottom boundary conditions" bottom_bc::BBC = LandComponentBC() "lateral boundary conditions" - lateral_bc::LBC = LandComponentBC() + minx_bc::LBC = LandComponentBC() + maxx_bc::LBC = LandComponentBC() + miny_bc::LBC = LandComponentBC() + maxy_bc::LBC = LandComponentBC() end """ @@ -128,9 +131,11 @@ needs to define the BC for the components they wish to model. Base.@kwdef struct LandComponentBC{ SW <: AbstractBoundaryConditions, SH <: AbstractBoundaryConditions, + R <: AbstractBoundaryConditions, } soil_water::SW = NoBC() soil_heat::SH = NoBC() + river::R = NoBC() end @@ -143,8 +148,8 @@ 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 diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 3792d36c67e..e5871dacb32 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -100,3 +100,106 @@ using ClimateMachine.BalanceLaws: @test size(Base.collect(keys(aux_vars)))[1] == 1 @test size(Base.collect(keys(state_vars)))[1] == 0 end + + + +@testset "Analytical River Model" begin + ClimateMachine.init() + FT = Float64 + + function init_land_model!(land, state, aux, localgeo, time) end + + 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_river = RiverModel( + slope_x = (x,y) -> eltype(x)(1), + slope_y = (x,y) -> eltype(x)(0.0), + mag_slope = (x,y) -> eltype(x)(0.0016), + width = (x,y) -> eltype(x)(1); + mannings = (x,y) -> eltype(x)(0.025) + ) + + sources = () + bc = LandDomainBC( + + ) + m = LandModel( + param_set, + m_soil, + m_river; + boundary_conditions = bc, + source = sources, + init_state_prognostic = init_land_model!, + ) + 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 + + N_poly = 1; + xres = FT(18.288) + 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.29) + driver_config = ClimateMachine.MultiColumnLandModel( + "LandModel", + (N_poly, N_poly), + (xres,yres,zres), + xmax, + ymax, + zmax, + param_set, + m; + zmin = zmin, + xmin = xmin, + ymin = ymin, + #numerical_flux_first_order = CentralNumericalFluxFirstOrder(),now the default for us + meshwarp = (x...) -> warp_constant_slope(x...; + topo_max = topo_max, zmin = zmin, xmax = xmax), + ); + + + t0 = FT(0) + timeend = FT(60) + 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 "z" as aux. + @test t == timeend + @test size(Base.collect(keys(aux_vars)))[1] == 1 + @test size(Base.collect(keys(state_vars)))[1] == 0 +end From 046aa36f97433f0a8b742591c8d935cd66b2b8db Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Wed, 13 Jan 2021 16:30:29 -0700 Subject: [PATCH 05/24] wip driver --- src/Land/Model/source.jl | 26 ++++++++++-- test/Land/Model/test_no_river.jl | 73 ++++++++++++++++++++------------ 2 files changed, 68 insertions(+), 31 deletions(-) diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index 8c004571f59..76463fd57b8 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -1,5 +1,5 @@ #### Land sources -export PhaseChange +export PhaseChange, Precip function heaviside(x::FT) where {FT} if x >= FT(0) @@ -11,14 +11,14 @@ function heaviside(x::FT) where {FT} end -abstract type SoilSource{FT <: AbstractFloat} end +abstract type LandSource{FT <: AbstractFloat} end """ PhaseChange <: SoilSource The function which computes the freeze/thaw source term for Richard's equation, assuming the timescale is the maximum of the thermal timescale and the timestep. """ -Base.@kwdef struct PhaseChange{FT} <: SoilSource{FT} +Base.@kwdef struct PhaseChange{FT} <: LandSource{FT} "Timestep" Δt::FT = FT(NaN) "Timescale for temperature changes" @@ -39,6 +39,26 @@ function land_source!( f(land, source, state, diffusive, aux, t, direction) end + +struct Precip{FT} <: LandSource{FT} end + precip::F +end +(p::Precip{FT})(x::FT, y::FT, t::FT) = FT(p.precip(x, y, t)) + + +function land_source!( + source_type::Precip, + land::LandModel, + source::Vars, + state::Vars, + diffusive::Vars, + aux::Vars, + t::Real, + direction, +) + source.river.height += source_type(aux.x, aux.y, t) +end + function land_source!( source_type::PhaseChange, land::LandModel, diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index e5871dacb32..8fbefe90f32 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -26,7 +26,7 @@ using ClimateMachine.VariableTemplates using ClimateMachine.SingleStackUtils using ClimateMachine.BalanceLaws: BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state - +#= @testset "NoRiver Model" begin ClimateMachine.init() FT = Float64 @@ -100,10 +100,9 @@ using ClimateMachine.BalanceLaws: @test size(Base.collect(keys(aux_vars)))[1] == 1 @test size(Base.collect(keys(state_vars)))[1] == 0 end +=# - - -@testset "Analytical River Model" begin +#@testset "Analytical River Model" begin ClimateMachine.init() FT = Float64 @@ -122,10 +121,22 @@ end mannings = (x,y) -> eltype(x)(0.025) ) - sources = () + + # 10.1061/(ASCE)0733-9429(2007)133:2(217) + # Eqn 6 bc = LandDomainBC( - + minx_bc = LandComponentBC(river = Dirichlet((aux,t) -> eltype(aux)(0))), ) + + function init_land_model!(land, state, aux, localgeo, time) + state.river.height = eltype(state)(0.0) + end + + # units in m / s + precip(x, y, t) = t < (30 * 60) ? 1.4e-5 : 0.0 + + sources = (Precip(precip),) + m = LandModel( param_set, m_soil, @@ -153,6 +164,7 @@ end xmax = FT(182.88) ymax = FT(1.0) topo_max = FT(0.29) + driver_config = ClimateMachine.MultiColumnLandModel( "LandModel", (N_poly, N_poly), @@ -170,9 +182,8 @@ end topo_max = topo_max, zmin = zmin, xmax = xmax), ); - t0 = FT(0) - timeend = FT(60) + timeend = FT(60*60) dt = FT(1) solver_config = ClimateMachine.SolverConfiguration( @@ -181,25 +192,31 @@ end driver_config, ode_dt = dt, ) - mygrid = solver_config.dg.grid + mygrid = solver_config.dg.grid Q = solver_config.Q - aux = solver_config.dg.state_auxiliary - - ClimateMachine.invoke!(solver_config;) + + height_index = + varsindex(vars_state(m, Prognostic(), FT), :river, :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) + height = Q[:, height_index, :] + all_vars = Dict{String, Array}( + "t" => [t], + "height" => height, + ) + dons[iostep[1]] = all_vars + iostep[1] += 1 + nothing + end - 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 "z" as aux. - @test t == timeend - @test size(Base.collect(keys(aux_vars)))[1] == 1 - @test size(Base.collect(keys(state_vars)))[1] == 0 -end + ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) +#end From 42dbda8af3532feee015cac9af9ea315b022127d Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Thu, 14 Jan 2021 14:06:04 -0700 Subject: [PATCH 06/24] wip river bc --- src/Land/Model/LandModel.jl | 4 ++ src/Land/Model/River.jl | 96 ++++++++++++++++++++++++-------- src/Land/Model/land_bc.jl | 13 +++-- src/Land/Model/source.jl | 14 +++-- test/Land/Model/test_no_river.jl | 39 +++++++------ 5 files changed, 116 insertions(+), 50 deletions(-) diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 3dbd724a5bc..8eebf58a937 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -96,6 +96,8 @@ end function vars_state(land::LandModel, st::Auxiliary, FT) @vars begin + x::FT + y::FT z::FT soil::vars_state(land.soil, st, FT) river::vars_state(land.river, st, FT) @@ -122,6 +124,8 @@ 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.river, aux, geom) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 0d83bd78e96..a450802c47e 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -1,24 +1,20 @@ module River -using ..VariableTemplates using DocStringExtensions using ..Land -using ..BalanceLaws +using ..VariableTemplates import ..BalanceLaws: BalanceLaw, vars_state, - flux_first_order!, - flux_second_order!, - source!, - boundary_conditions, - boundary_state!, - compute_gradient_argument!, - compute_gradient_flux!, - nodal_init_state_auxiliary!, - init_state_prognostic!, - nodal_update_auxiliary_state! - -export RiverModel, NoRiverModel + flux_first_order!, + Prognostic, + Auxiliary, + Gradient, + GradientFlux +using ...DGMethods: LocalGeometry +using StaticArrays: SVector + +export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state! struct NoRiverModel <: BalanceLaw end @@ -30,7 +26,7 @@ struct RiverModel{M,Sx,Sy,MS,W} <: BalanceLaw width::W end -function axyerModel( +function RiverModel( slope_x::Function, slope_y::Function, mag_slope::Function, @@ -54,12 +50,12 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) return SVector(sx * magnitude, sy * magnitude, zero(FT)) end -vars_state(water::RiverModel, st::Prognostic, FT) = @vars(height::FT) -vars_state(water::RiverModel, st::Auxiliary, FT) = @vars() -vars_state(water::RiverModel, st::Gradient, FT) = @vars() -vars_state(water::RiverModel, st::GradientFlux, FT) = @vars() +vars_state(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) +vars_state(river::RiverModel, st::Auxiliary, FT) = @vars() +vars_state(river::RiverModel, st::Gradient, FT) = @vars() +vars_state(river::RiverModel, st::GradientFlux, FT) = @vars() -function Land.land_init_aux!(land::LandModel, river::BalanceLaw, aux, geom) +function Land.land_init_aux!(land::LandModel, river::BalanceLaw, aux, geom::LocalGeometry) end function Land.compute_gradient_argument!(land::LandModel, river::BalanceLaw, transform::Grad, state, aux, t) @@ -73,13 +69,67 @@ end function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) - x, y = aux.x, aux.y + x = aux.x + y = aux.y width = river.width(x, y) height = state.river.area / width - v = calculate_velocity(river, x, y ,height) + v = calculate_velocity(river, x, y, height) Q = state.river.area * v - flux.river.height = Q + flux.river.area = Q end +# boundry conditions + +# General case - to be used with bc::NoBC +function river_boundary_flux!( + nf, + bc::Land.AbstractBoundaryConditions, + m, + land::LandModel, + _..., +) +end + +function river_boundary_state!( + nf, + bc::Land.AbstractBoundaryConditions, + m, + land::LandModel, + _..., +) +end + +# Dirichlet BC for River +function river_boundary_flux!( + nf, + bc::Land.Dirichlet, + model::RiverModel, + land::LandModel, + state⁺::Vars, + aux⁺::Vars, + nM, + state⁻::Vars, + aux⁻::Vars, + t, + _..., +) +end + +function river_boundary_state!( + nf, + bc::Land.Dirichlet, + model::RiverModel, + land::LandModel, + state⁺::Vars, + aux⁺::Vars, + nM, + state⁻::Vars, + aux⁻::Vars, + t, + _..., +) + bc_function = bc.state_bc + state⁺.river.area = bc_function(aux⁻, t) +end end \ No newline at end of file diff --git a/src/Land/Model/land_bc.jl b/src/Land/Model/land_bc.jl index e3f68f08f82..d1699a78dc6 100644 --- a/src/Land/Model/land_bc.jl +++ b/src/Land/Model/land_bc.jl @@ -106,16 +106,16 @@ 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" - minx_bc::LBC = LandComponentBC() - maxx_bc::LBC = LandComponentBC() - miny_bc::LBC = LandComponentBC() - maxy_bc::LBC = LandComponentBC() + minx_bc::MinXBC = LandComponentBC() + maxx_bc::MaxXBC = LandComponentBC() + miny_bc::MinYBC = LandComponentBC() + maxy_bc::MaxYBC = LandComponentBC() end """ @@ -181,10 +181,10 @@ function boundary_state!( ) 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...) + river_boundary_state!(nf, bc.river, land.river, land, args...) end @@ -222,4 +222,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...) + river_boundary_flux!(nf, bc.river, land.river, land, args...) end diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index 76463fd57b8..5c80c079c1e 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -39,12 +39,18 @@ function land_source!( f(land, source, state, diffusive, aux, t, direction) end - -struct Precip{FT} <: LandSource{FT} end +#TODO: precip docs +struct Precip{FT, F} <: LandSource{FT} precip::F + + function Precip{FT}(precip::F) where {FT, F} + new{FT, F}(precip) + end end -(p::Precip{FT})(x::FT, y::FT, t::FT) = FT(p.precip(x, y, t)) +function (p::Precip{FT})(x, y, t) where {FT} + FT(p.precip(x, y, t)) +end function land_source!( source_type::Precip, @@ -56,7 +62,7 @@ function land_source!( t::Real, direction, ) - source.river.height += source_type(aux.x, aux.y, t) + source.river.area += source_type(aux.x, aux.y, t) end function land_source!( diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 8fbefe90f32..0d99ce78d6c 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -114,28 +114,29 @@ end m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) m_river = RiverModel( - slope_x = (x,y) -> eltype(x)(1), - slope_y = (x,y) -> eltype(x)(0.0), - mag_slope = (x,y) -> eltype(x)(0.0016), - width = (x,y) -> eltype(x)(1); + (x,y) -> eltype(x)(1), + (x,y) -> eltype(x)(0.0), + (x,y) -> eltype(x)(0.0016), + (x,y) -> eltype(x)(1); mannings = (x,y) -> eltype(x)(0.025) ) - # 10.1061/(ASCE)0733-9429(2007)133:2(217) # Eqn 6 + # debug boundry condition, constant positive flow rate after min 30 + # add river boundry state land domain bc bc = LandDomainBC( - minx_bc = LandComponentBC(river = Dirichlet((aux,t) -> eltype(aux)(0))), + minx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), ) function init_land_model!(land, state, aux, localgeo, time) - state.river.height = eltype(state)(0.0) + state.river.area = eltype(state)(0) end # units in m / s precip(x, y, t) = t < (30 * 60) ? 1.4e-5 : 0.0 - sources = (Precip(precip),) + sources = (Precip{FT}(precip),) m = LandModel( param_set, @@ -175,8 +176,6 @@ end param_set, m; zmin = zmin, - xmin = xmin, - ymin = ymin, #numerical_flux_first_order = CentralNumericalFluxFirstOrder(),now the default for us meshwarp = (x...) -> warp_constant_slope(x...; topo_max = topo_max, zmin = zmin, xmax = xmax), @@ -192,11 +191,11 @@ end driver_config, ode_dt = dt, ) - mygrid = solver_config.dg.grid + mygrid = solver_config.dg.grid Q = solver_config.Q - - height_index = - varsindex(vars_state(m, Prognostic(), FT), :river, :height) + + area_index = + varsindex(vars_state(m, Prognostic(), FT), :river, :area) n_outputs = 60 every_x_simulation_time = ceil(Int, timeend / n_outputs) @@ -208,15 +207,21 @@ end every_x_simulation_time, ) do (init = false) t = ODESolvers.gettime(solver_config.solver) - height = Q[:, height_index, :] + area = Q[:, area_index, :] all_vars = Dict{String, Array}( "t" => [t], - "height" => height, + "area" => area, ) dons[iostep[1]] = all_vars iostep[1] += 1 - nothing + return end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) + +# using Statistics +# mask = aux[:,-1,:] .== 182.88 +# area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] +# velocity = area.^(5/3) .* (sqrt(0.0016)/0.025) + #end From d4109fa179198b9098155f31e64110745ceb8bdf Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Thu, 14 Jan 2021 15:38:23 -0700 Subject: [PATCH 07/24] wip river bc --- src/Land/Model/River.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index a450802c47e..5d078fd2673 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -18,12 +18,12 @@ export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state! struct NoRiverModel <: BalanceLaw end -struct RiverModel{M,Sx,Sy,MS,W} <: BalanceLaw - mannings::M +struct RiverModel{Sx,Sy,MS,W,M} <: BalanceLaw slope_x::Sx slope_y::Sy mag_slope::MS width::W + mannings::M end function RiverModel( From d137568733e1b6472aa61e4c2bb7ee8f52696c37 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 14 Jan 2021 15:51:24 -0800 Subject: [PATCH 08/24] added analytic function to test --- test/Land/Model/test_no_river.jl | 54 +++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 0d99ce78d6c..6e339c008dd 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -5,6 +5,7 @@ using MPI using OrderedCollections using StaticArrays using Test +using Statistics using CLIMAParameters struct EarthParameterSet <: AbstractEarthParameterSet end @@ -219,9 +220,54 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) -# using Statistics -# mask = aux[:,-1,:] .== 182.88 -# area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] -# velocity = area.^(5/3) .* (sqrt(0.0016)/0.025) +aux = solver_config.dg.state_auxiliary; +mask = aux[:,1,:] .== 182.88 +area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] +velocity = area.^(5/3) .* (sqrt(0.0016)/0.025) #m /s (this isnt right??) +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 + + +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 + + #end From 8568e64fea4cf4a8c94c040e5fab6bfd03e3f596 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Fri, 15 Jan 2021 08:03:27 -0800 Subject: [PATCH 09/24] both no river and analytic case 1 run --- src/Land/Model/River.jl | 21 ++++--- test/Land/Model/test_no_river.jl | 97 ++++++++++++++++---------------- 2 files changed, 60 insertions(+), 58 deletions(-) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 5d078fd2673..735740afb03 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -11,11 +11,13 @@ import ..BalanceLaws: Auxiliary, Gradient, GradientFlux + using ...DGMethods: LocalGeometry using StaticArrays: SVector export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state! + struct NoRiverModel <: BalanceLaw end struct RiverModel{Sx,Sy,MS,W,M} <: BalanceLaw @@ -46,29 +48,26 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) FT = eltype(h) sx = FT(river.slope_x(x, y)) sy = FT(river.slope_y(x, y)) - magnitude = h^FT(2/3) / (river.mannings(x, y) * sqrt(river.mag_slope(x, y))) + magnitude = h^FT(2/3) / river.mannings(x, y) * sqrt(river.mag_slope(x, y)) return SVector(sx * magnitude, sy * magnitude, zero(FT)) end +## there is a default method for balance laws that adds no variables that we can use, so +### i deleted the other methods. vars_state(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) -vars_state(river::RiverModel, st::Auxiliary, FT) = @vars() -vars_state(river::RiverModel, st::Gradient, FT) = @vars() -vars_state(river::RiverModel, st::GradientFlux, FT) = @vars() function Land.land_init_aux!(land::LandModel, river::BalanceLaw, aux, geom::LocalGeometry) end -function Land.compute_gradient_argument!(land::LandModel, river::BalanceLaw, transform::Grad, state, aux, t) - # v = calculate_velocity(river, state.river.height,river) - # transform.river.grad_Q = state.river.height*v#*river.width(aux.x,aux.y) -end - function Land.land_nodal_update_auxiliary_state!(land::LandModel, river::BalanceLaw, state, aux, t) end +function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +end + -function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) x = aux.x y = aux.y width = river.width(x, y) @@ -132,4 +131,4 @@ function river_boundary_state!( state⁺.river.area = bc_function(aux⁻, t) end -end \ No newline at end of file +end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 6e339c008dd..7dbb08dd5bc 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -27,7 +27,7 @@ using ClimateMachine.VariableTemplates using ClimateMachine.SingleStackUtils using ClimateMachine.BalanceLaws: BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state -#= + @testset "NoRiver Model" begin ClimateMachine.init() FT = Float64 @@ -98,12 +98,12 @@ using ClimateMachine.BalanceLaws: ) #Make sure it runs, and that there are no state variables, and only "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 -=# -#@testset "Analytical River Model" begin + +@testset "Analytical River Model" begin ClimateMachine.init() FT = Float64 @@ -157,7 +157,7 @@ end end N_poly = 1; - xres = FT(18.288) + xres = FT(2.286) yres = FT(0.25) zres = FT(0.1) # Specify the domain boundaries. @@ -165,7 +165,7 @@ end zmin = FT(-0.1); xmax = FT(182.88) ymax = FT(1.0) - topo_max = FT(0.29) + topo_max = FT(0.0016*xmax) driver_config = ClimateMachine.MultiColumnLandModel( "LandModel", @@ -221,53 +221,56 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) -aux = solver_config.dg.state_auxiliary; -mask = aux[:,1,:] .== 182.88 -area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] -velocity = area.^(5/3) .* (sqrt(0.0016)/0.025) #m /s (this isnt right??) -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 - - -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) + aux = solver_config.dg.state_auxiliary; + mask = aux[:,1,:] .== 182.88 + n_outputs = length(dons) + area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] + height = area ./ ymax + 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 - if t <= t_r && t > t_c - return alpha*(i*t_c)^(m) + 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 - - 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 + + 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 - return alpha*yL^m end + @test sqrt_rmse_over_max_q = sqrt(mean((analytic.(time_data, Ref(alpha), Ref(t_c), Ref(t_r), Ref(i), Ref(L), Ref(m)) .- q).^2.0))/ maximum(q) < 3e-3 + end - - -#end From 25501ae28ded87c9d1773973b61dbaa38e9edd16 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Thu, 21 Jan 2021 10:55:52 -0700 Subject: [PATCH 10/24] wip cleanup --- src/Land/Model/River.jl | 14 ++++------ test/Land/Model/test_no_river.jl | 47 ++++++++++++++++---------------- 2 files changed, 30 insertions(+), 31 deletions(-) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 735740afb03..4db20cb9e4c 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -56,18 +56,16 @@ end ### i deleted the other methods. vars_state(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) -function Land.land_init_aux!(land::LandModel, river::BalanceLaw, aux, geom::LocalGeometry) +function Land.land_init_aux!(land::LandModel, river::Union{NoRiverModel,RiverModel}, aux, geom::LocalGeometry) end - -function Land.land_nodal_update_auxiliary_state!(land::LandModel, river::BalanceLaw, state, aux, t) +function Land.land_nodal_update_auxiliary_state!(land::LandModel, river::Union{NoRiverModel,RiverModel}, state, aux, t) end -function flux_first_order!(land::LandModel, river::BalanceLaw, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +function flux_first_order!(land::LandModel, river::NoRiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) end - -function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) x = aux.x y = aux.y width = river.width(x, y) @@ -83,7 +81,7 @@ end function river_boundary_flux!( nf, bc::Land.AbstractBoundaryConditions, - m, + m::Union{NoRiverModel,RiverModel}, land::LandModel, _..., ) @@ -92,7 +90,7 @@ end function river_boundary_state!( nf, bc::Land.AbstractBoundaryConditions, - m, + m::Union{NoRiverModel,RiverModel}, land::LandModel, _..., ) diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 7dbb08dd5bc..614ec6a77ed 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -28,12 +28,13 @@ using ClimateMachine.SingleStackUtils using ClimateMachine.BalanceLaws: BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state +function init_land_model!(land, state, aux, localgeo, time) +end + @testset "NoRiver Model" begin ClimateMachine.init() FT = Float64 - function init_land_model!(land, state, aux, localgeo, time) end - soil_water_model = PrescribedWaterModel() soil_heat_model = PrescribedTemperatureModel() soil_param_functions = nothing @@ -68,7 +69,6 @@ using ClimateMachine.BalanceLaws: numerical_flux_first_order = CentralNumericalFluxFirstOrder(), ) - t0 = FT(0) timeend = FT(60) dt = FT(1) @@ -102,13 +102,19 @@ using ClimateMachine.BalanceLaws: @test size(Base.collect(keys(state_vars)))[1] == 0 end +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 @testset "Analytical River Model" begin ClimateMachine.init() FT = Float64 - function init_land_model!(land, state, aux, localgeo, time) end - soil_water_model = PrescribedWaterModel() soil_heat_model = PrescribedTemperatureModel() soil_param_functions = nothing @@ -121,11 +127,10 @@ end (x,y) -> eltype(x)(1); mannings = (x,y) -> eltype(x)(0.025) ) - + + # Analytical test case defined as Model 1 in DOI: # 10.1061/(ASCE)0733-9429(2007)133:2(217) # Eqn 6 - # debug boundry condition, constant positive flow rate after min 30 - # add river boundry state land domain bc bc = LandDomainBC( minx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), ) @@ -147,14 +152,6 @@ end source = sources, init_state_prognostic = init_land_model!, ) - 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 N_poly = 1; xres = FT(2.286) @@ -177,7 +174,6 @@ end param_set, m; zmin = zmin, - #numerical_flux_first_order = CentralNumericalFluxFirstOrder(),now the default for us meshwarp = (x...) -> warp_constant_slope(x...; topo_max = topo_max, zmin = zmin, xmax = xmax), ); @@ -220,14 +216,19 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) - + # Compare flowrate analytical derivation aux = solver_config.dg.state_auxiliary; - mask = aux[:,1,:] .== 182.88 + + # get all nodal points at the max X bound of the domain + mask = Array(aux[:,1,:] .== 182.88) n_outputs = length(dons) - area = [mean(dons[k]["area"][mask[:]]) for k in 1:n_outputs] + # get prognostic variable area from nodal state (m^2) + area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] height = area ./ ymax + # 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 @@ -271,6 +272,6 @@ end end - @test sqrt_rmse_over_max_q = sqrt(mean((analytic.(time_data, Ref(alpha), Ref(t_c), Ref(t_r), Ref(i), Ref(L), Ref(m)) .- q).^2.0))/ maximum(q) < 3e-3 - + # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a GPU array) + @test sqrt_rmse_over_max_q = sqsrt(mean((analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q).^4.0))/ maximum(q) < 3e-3 end From 3824fdf33ed8657f14745b980e928141ee8556f2 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Thu, 21 Jan 2021 13:14:12 -0700 Subject: [PATCH 11/24] wip maxwell v catchment --- src/Land/Model/River.jl | 18 ++-- test/Land/Model/test_no_river.jl | 163 ++++++++++++++++++++++++++++++- 2 files changed, 171 insertions(+), 10 deletions(-) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 4db20cb9e4c..e3c36f3ff0f 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -19,11 +19,11 @@ export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state! struct NoRiverModel <: BalanceLaw end - -struct RiverModel{Sx,Sy,MS,W,M} <: BalanceLaw +## we should change this - we should only specify the slope vector components, not the unit vector components +## and the magnitude +struct RiverModel{Sx,Sy,W,M} <: BalanceLaw slope_x::Sx slope_y::Sy - mag_slope::MS width::W mannings::M end @@ -31,15 +31,13 @@ end function RiverModel( slope_x::Function, slope_y::Function, - mag_slope::Function, width::Function; mannings::Function = (x, y) -> convert(eltype(x), 0.03)) args = ( slope_x, slope_y, - mag_slope, width, - mannings + mannings, ) return RiverModel{typeof.(args)...}(args...) end @@ -48,8 +46,12 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) FT = eltype(h) sx = FT(river.slope_x(x, y)) sy = FT(river.slope_y(x, y)) - magnitude = h^FT(2/3) / river.mannings(x, y) * sqrt(river.mag_slope(x, y)) - return SVector(sx * magnitude, sy * magnitude, zero(FT)) + magnitude_slope = FT(sqrt(sx^FT(2.0)+sy^FT(2.0))) + mannings_coeff = FT(river.mannings(x, y)) + magnitude = h^FT(2/3) / mannings_coeff * sqrt(magnitude_slope) + return SVector(sx * magnitude / magnitude_slope, + sy * magnitude / magnitude_slope, + zero(FT)) end ## there is a default method for balance laws that adds no variables that we can use, so diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 614ec6a77ed..6f17daaf291 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -31,6 +31,7 @@ using ClimateMachine.BalanceLaws: function init_land_model!(land, state, aux, localgeo, time) end +#= @testset "NoRiver Model" begin ClimateMachine.init() FT = Float64 @@ -111,6 +112,7 @@ function warp_constant_slope(xin, yin, zin; topo_max = 0.2, zmin = -0.1, xmax = return x, y, z end + @testset "Analytical River Model" begin ClimateMachine.init() FT = Float64 @@ -120,10 +122,11 @@ end soil_param_functions = nothing m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) + m_river = RiverModel( - (x,y) -> eltype(x)(1), - (x,y) -> eltype(x)(0.0), (x,y) -> eltype(x)(0.0016), + (x,y) -> eltype(x)(0.0), +# (x,y) -> eltype(x)(0.0016), (x,y) -> eltype(x)(1); mannings = (x,y) -> eltype(x)(0.025) ) @@ -275,3 +278,159 @@ end # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a GPU array) @test sqrt_rmse_over_max_q = sqsrt(mean((analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q).^4.0))/ maximum(q) < 3e-3 end +=# + +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 + +## @testset "V Catchment Maxwell River Model" begin + 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 * 10^-4) : MFT(2.5 * 10^-3) + end + + m_river = RiverModel( + x_slope, + y_slope, + (x,y) -> eltype(x)(1); + mannings = channel_mannings, + ) + + bc = LandDomainBC( + miny_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), + ) + + function init_land_model!(land, state, aux, localgeo, time) + state.river.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, + m_river; + boundary_conditions = bc, + source = sources, + init_state_prognostic = init_land_model!, + ) + + N_poly = 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, N_poly), + (xres,yres,zres), + xmax, + ymax, + zmax, + param_set, + m; + zmin = zmin, + meshwarp = (x...) -> warp_tilted_v(x...), + ); + + t0 = FT(0) + timeend = FT(180*60) + dt = FT(30) + + 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), :river, :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,)) + + # 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 area from nodal state (m^2) + area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] + height = area ./ ymax + + # get similation timesteps (s) + time_data = [dons[l]["t"][1] for l in 1:n_outputs] +#end \ No newline at end of file From 214e071f4db0fedc0e0eda31097c4f71cfa2c862 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 28 Jan 2021 16:40:45 -0800 Subject: [PATCH 12/24] still wip maxwell tilted V --- src/Land/Model/River.jl | 18 +++--- test/Land/Model/test_no_river.jl | 94 +++++++++++++++++++------------- 2 files changed, 64 insertions(+), 48 deletions(-) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index e3c36f3ff0f..04469745f3f 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -15,7 +15,7 @@ import ..BalanceLaws: using ...DGMethods: LocalGeometry using StaticArrays: SVector -export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state! +export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state!, calculate_velocity struct NoRiverModel <: BalanceLaw end @@ -46,16 +46,15 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) FT = eltype(h) sx = FT(river.slope_x(x, y)) sy = FT(river.slope_y(x, y)) - magnitude_slope = FT(sqrt(sx^FT(2.0)+sy^FT(2.0))) + #magnitude_slope = FT(sqrt(sx^FT(2.0)+sy^FT(2.0))) mannings_coeff = FT(river.mannings(x, y)) - magnitude = h^FT(2/3) / mannings_coeff * sqrt(magnitude_slope) - return SVector(sx * magnitude / magnitude_slope, - sy * magnitude / magnitude_slope, + magnitude = h^FT(2/3) / mannings_coeff #* sqrt(magnitude_slope) + #if the slope is positive, dz/dx >0, flow should be in opposite direction. add in minus signs + return SVector(-sign(sx) * magnitude*sqrt(abs(sx)),# / magnitude_slope, + -sign(sy) * magnitude*sqrt(abs(sy)),# / magnitude_slope, zero(FT)) end -## there is a default method for balance laws that adds no variables that we can use, so -### i deleted the other methods. vars_state(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) function Land.land_init_aux!(land::LandModel, river::Union{NoRiverModel,RiverModel}, aux, geom::LocalGeometry) @@ -71,9 +70,10 @@ function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state x = aux.x y = aux.y width = river.width(x, y) - height = state.river.area / width + area = max(eltype(state)(0.0), state.river.area) + height = area / width v = calculate_velocity(river, x, y, height) - Q = state.river.area * v + Q = area * v flux.river.area = Q end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 6f17daaf291..cace7316891 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -31,7 +31,7 @@ using ClimateMachine.BalanceLaws: function init_land_model!(land, state, aux, localgeo, time) end -#= + @testset "NoRiver Model" begin ClimateMachine.init() FT = Float64 @@ -124,7 +124,7 @@ end m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) m_river = RiverModel( - (x,y) -> eltype(x)(0.0016), + (x,y) -> eltype(x)(-0.0016), (x,y) -> eltype(x)(0.0), # (x,y) -> eltype(x)(0.0016), (x,y) -> eltype(x)(1); @@ -276,9 +276,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) - @test sqrt_rmse_over_max_q = sqsrt(mean((analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q).^4.0))/ maximum(q) < 3e-3 + @test sqrt_rmse_over_max_q = sqrt(mean((analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q).^4.0))/ maximum(q) < 3e-3 end -=# + function warp_tilted_v(xin, yin, zin) FT = eltype(xin) @@ -298,19 +298,19 @@ function warp_tilted_v(xin, yin, zin) return x, y, z end -## @testset "V Catchment Maxwell River Model" begin +#@testset "V Catchment Maxwell River Model" begin 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) + if x < MFT(800) MFT(-0.05) elseif x <= MFT(820) MFT(0) @@ -318,37 +318,39 @@ end 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 * 10^-4) : MFT(2.5 * 10^-3) + return x >= MFT(800) && x <= MFT(820) ? MFT(2.5*60 * 10^-3) : MFT(2.5*60 * 10^-4) end - + m_river = RiverModel( x_slope, y_slope, (x,y) -> eltype(x)(1); mannings = channel_mannings, ) - + bc = LandDomainBC( miny_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), + minx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), + maxx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), ) - + function init_land_model!(land, state, aux, localgeo, time) state.river.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, @@ -357,8 +359,9 @@ end source = sources, init_state_prognostic = init_land_model!, ) - - N_poly = 1; + + N_poly_hori = 1; + N_poly_vert = 1; xres = FT(20) yres = FT(20) zres = FT(1) @@ -367,11 +370,11 @@ end zmin = FT(0); xmax = FT(1620) ymax = FT(1000) - - + + driver_config = ClimateMachine.MultiColumnLandModel( "LandModel", - (N_poly, N_poly), + (N_poly_hori, N_poly_vert), (xres,yres,zres), xmax, ymax, @@ -379,13 +382,13 @@ end param_set, m; zmin = zmin, - meshwarp = (x...) -> warp_tilted_v(x...), + # meshwarp = (x...) -> warp_tilted_v(x...), ); - + t0 = FT(0) timeend = FT(180*60) - dt = FT(30) - + dt = FT(3) + solver_config = ClimateMachine.SolverConfiguration( t0, timeend, @@ -398,11 +401,11 @@ end area_index = varsindex(vars_state(m, Prognostic(), FT), :river, :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, @@ -417,20 +420,33 @@ end iostep[1] += 1 return end - + ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) - + # Compare flowrate analytical derivation aux = solver_config.dg.state_auxiliary; - + x = aux[:,1,:] + y = aux[:,2,:] + z = aux[:,3,:] # 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 area from nodal state (m^2) - area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] - height = area ./ ymax + mask2 = (Float64.(aux[:,2,:] .== 1000.0) .+ Float64.(aux[:,1,:] .<= 820) .+ Float64.(aux[:,1,:] .> 800)) .== 3 + n_outputs = length(dons) + function compute_Q(a,xv) + height = max.(a,0.0) ./ 1.0# width = 1 here, everywhere + v = calculate_velocity(m_river,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 + # Maybe there is a conceptual error here? + Q = [sum(compute_Q.(Array(dons[k]["area"])[:][mask2[:]], x[mask2[:]])) for k in 1:n_outputs] ./4.0 .* xres + + # get similation timesteps (s) time_data = [dons[l]["t"][1] for l in 1:n_outputs] -#end \ No newline at end of file + #also helpful + #scatter(aux[:,1,:][:], aux[:,2,:][:], dons[50]["area"][:]) +#end From 77a7e822ce5ddb2a8d87c3ce2b6ecfbb91d1aa64 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Wed, 3 Feb 2021 17:12:35 -0800 Subject: [PATCH 13/24] added in wavespeed method for Rusanov flux --- src/Land/Model/LandModel.jl | 25 ++++++++++++++++++++++++- src/Land/Model/River.jl | 1 + test/Land/Model/test_no_river.jl | 7 ++++--- 3 files changed, 29 insertions(+), 4 deletions(-) diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 8eebf58a937..a58baebc5af 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -22,7 +22,8 @@ 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 @@ -233,6 +234,9 @@ function init_state_prognostic!( land.init_state_prognostic(land, state, aux, coords, t, args...) end + + + include("Runoff.jl") using .Runoff include("land_bc.jl") @@ -247,4 +251,23 @@ include("soil_bc.jl") include("source.jl") include("River.jl") using .River + +function wavespeed( + m::LandModel, + n⁻, + state::Vars, + aux::Vars, + t::Real, + direction +) + FT = eltype(state) + g = FT(9.8) + width = m.river.width(aux.x,aux.y) + area = max(eltype(state)(0.0), state.river.area) + height = area ./ width + v = calculate_velocity(m.river, aux.x, aux.y, height) + speed = abs(norm(v)) + return speed#+sqrt(g*height) +end + end # Module diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 04469745f3f..767a0fdb052 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -131,4 +131,5 @@ function river_boundary_state!( state⁺.river.area = bc_function(aux⁻, t) end + end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index cace7316891..87634058a04 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -183,7 +183,7 @@ end t0 = FT(0) timeend = FT(60*60) - dt = FT(1) + dt = FT(10) solver_config = ClimateMachine.SolverConfiguration( t0, @@ -382,12 +382,13 @@ end 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(3) + dt = FT(0.5) solver_config = ClimateMachine.SolverConfiguration( t0, @@ -429,7 +430,7 @@ end y = aux[:,2,:] z = aux[:,3,:] # get all nodal points at the max X bound of the domain - mask2 = (Float64.(aux[:,2,:] .== 1000.0) .+ Float64.(aux[:,1,:] .<= 820) .+ Float64.(aux[:,1,:] .> 800)) .== 3 + mask2 = (Float64.(aux[:,2,:] .== 1000.0)) .==1# .+ Float64.(aux[:,1,:] .<= 820) .+ Float64.(aux[:,1,:] .> 800)) .== 3 n_outputs = length(dons) function compute_Q(a,xv) height = max.(a,0.0) ./ 1.0# width = 1 here, everywhere From aaee95661bf12b7d602457c7c264687b9cf1e58c Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 10:45:37 -0700 Subject: [PATCH 14/24] changed names, documentation, added tilted V artifact, kept old River stuff --- docs/list_of_apis.jl | 1 + docs/src/APIs/Land/SurfaceFlow.md | 14 + src/Land/Model/LandModel.jl | 38 +-- src/Land/Model/River.jl | 14 +- src/Land/Model/SurfaceFlow.jl | 182 ++++++++++ src/Land/Model/land_bc.jl | 10 +- src/Land/Model/source.jl | 2 +- test/Land/Model/Artifacts.toml | 3 + test/Land/Model/runtests.jl | 2 +- test/Land/Model/test_no_river.jl | 143 ++++---- test/Land/Model/test_overland_flow.jl | 456 ++++++++++++++++++++++++++ 11 files changed, 760 insertions(+), 105 deletions(-) create mode 100644 docs/src/APIs/Land/SurfaceFlow.md create mode 100644 src/Land/Model/SurfaceFlow.jl create mode 100644 test/Land/Model/test_overland_flow.jl diff --git a/docs/list_of_apis.jl b/docs/list_of_apis.jl index a9e835af935..989bbf0b703 100644 --- a/docs/list_of_apis.jl +++ b/docs/list_of_apis.jl @@ -22,6 +22,7 @@ apis = [ "APIs/Land/SoilWaterParameterizations.md", "Soil Heat Parameterizations" => "APIs/Land/SoilHeatParameterizations.md", + "Surface Flow" => "APIs/Land/SurfaceFlow.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..ac8beca7920 --- /dev/null +++ b/docs/src/APIs/Land/SurfaceFlow.md @@ -0,0 +1,14 @@ +# SurfaceFlow + +```@meta +CurrentModule = ClimateMachine.Land.SurfaceFlow +``` + +## SurfaceFlow types and functions +```@docs +OverlandFlowModel +NoSurfaceFlowModel +surface_boundary_flux! +surface_boundary_state! +calculate_velocity +``` \ No newline at end of file diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index a58baebc5af..66f9580d7b1 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -47,13 +47,13 @@ Users may over-ride prescribed default values for each field. # Fields $(DocStringExtensions.FIELDS) """ -struct LandModel{PS, S, R, LBC, SRC, IS} <: BalanceLaw +struct LandModel{PS, S, SF, LBC, SRC, IS} <: BalanceLaw "Parameter set" param_set::PS "Soil model" soil::S - "River model" - river::R + "Surface Flow model" + surface::SF "struct of boundary conditions" boundary_conditions::LBC "Source Terms (Problem specific source terms)" @@ -76,13 +76,13 @@ Constructor for the LandModel structure. function LandModel( param_set::AbstractParameterSet, soil::BalanceLaw, - river::BalanceLaw; + surface::BalanceLaw; boundary_conditions::LBC = LandDomainBC(), source::SRC = (), init_state_prognostic::IS = nothing, ) where {SRC, IS, LBC} @assert init_state_prognostic ≠ nothing - land = (param_set, soil, river, boundary_conditions, source, init_state_prognostic) + land = (param_set, soil, surface, boundary_conditions, source, init_state_prognostic) return LandModel{typeof.(land)...}(land...) end @@ -90,7 +90,7 @@ end function vars_state(land::LandModel, st::Prognostic, FT) @vars begin soil::vars_state(land.soil, st, FT) - river::vars_state(land.river, st, FT) + surface::vars_state(land.surface, st, FT) end end @@ -101,21 +101,21 @@ function vars_state(land::LandModel, st::Auxiliary, FT) y::FT z::FT soil::vars_state(land.soil, st, FT) - river::vars_state(land.river, 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) - river::vars_state(land.river, 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) - river::vars_state(land.river, st, FT) + surface::vars_state(land.surface, st, FT) end end @@ -129,7 +129,7 @@ function nodal_init_state_auxiliary!( aux.y = geom.coord[2] aux.z = geom.coord[3] land_init_aux!(land, land.soil, aux, geom) - land_init_aux!(land, land.river, aux, geom) + land_init_aux!(land, land.surface, aux, geom) end function flux_first_order!( @@ -140,7 +140,7 @@ function flux_first_order!( t::Real, directions, ) - flux_first_order!(land, land.river, flux, state, aux, t, directions) + flux_first_order!(land, land.surface, flux, state, aux, t, directions) end @@ -153,7 +153,7 @@ function compute_gradient_argument!( ) compute_gradient_argument!(land, land.soil, transform, state, aux, t) - compute_gradient_argument!(land, land.river, transform, state, aux, t) + compute_gradient_argument!(land, land.surface, transform, state, aux, t) end function compute_gradient_flux!( @@ -206,7 +206,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.river, state, aux, t) + land_nodal_update_auxiliary_state!(land, land.surface, state, aux, t) end @@ -249,8 +249,8 @@ include("soil_water.jl") include("soil_heat.jl") include("soil_bc.jl") include("source.jl") -include("River.jl") -using .River +include("SurfaceFlow.jl") +using .SurfaceFlow function wavespeed( m::LandModel, @@ -262,12 +262,12 @@ function wavespeed( ) FT = eltype(state) g = FT(9.8) - width = m.river.width(aux.x,aux.y) - area = max(eltype(state)(0.0), state.river.area) + width = m.surface.width(aux.x,aux.y) + area = max(eltype(state)(0.0), state.surface.area) height = area ./ width - v = calculate_velocity(m.river, aux.x, aux.y, height) + v = calculate_velocity(m.surface, aux.x, aux.y, height) speed = abs(norm(v)) - return speed#+sqrt(g*height) + return speed end end # Module diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index 767a0fdb052..c7058c722eb 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -17,10 +17,7 @@ using StaticArrays: SVector export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state!, calculate_velocity - struct NoRiverModel <: BalanceLaw end -## we should change this - we should only specify the slope vector components, not the unit vector components -## and the magnitude struct RiverModel{Sx,Sy,W,M} <: BalanceLaw slope_x::Sx slope_y::Sy @@ -46,12 +43,11 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) FT = eltype(h) sx = FT(river.slope_x(x, y)) sy = FT(river.slope_y(x, y)) - #magnitude_slope = FT(sqrt(sx^FT(2.0)+sy^FT(2.0))) mannings_coeff = FT(river.mannings(x, y)) - magnitude = h^FT(2/3) / mannings_coeff #* sqrt(magnitude_slope) - #if the slope is positive, dz/dx >0, flow should be in opposite direction. add in minus signs - return SVector(-sign(sx) * magnitude*sqrt(abs(sx)),# / magnitude_slope, - -sign(sy) * magnitude*sqrt(abs(sy)),# / magnitude_slope, + 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 @@ -77,7 +73,7 @@ function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state flux.river.area = Q end -# boundry conditions +# Boundary Conditions # General case - to be used with bc::NoBC function river_boundary_flux!( diff --git a/src/Land/Model/SurfaceFlow.jl b/src/Land/Model/SurfaceFlow.jl new file mode 100644 index 00000000000..1b9f6a119bf --- /dev/null +++ b/src/Land/Model/SurfaceFlow.jl @@ -0,0 +1,182 @@ +module SurfaceFlow + +using DocStringExtensions +using ..Land +using ..VariableTemplates +import ..BalanceLaws: + BalanceLaw, + vars_state, + flux_first_order!, + Prognostic, + Auxiliary, + Gradient, + GradientFlux + +using ...DGMethods: LocalGeometry +using StaticArrays: SVector + +export OverlandFlowModel, NoSurfaceFlowModel, surface_boundary_flux!, surface_boundary_state!, calculate_velocity + +""" + 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,W,M} <: BalanceLaw + +The 2D overland flow model - this model adds a prognostic variable +to the LandModel reflecting the surface area of water per unit length. +This model simulates depth averaged shallow water equation under the +kinematic approximation, along with Manning's relationship. +""" +struct OverlandFlowModel{Sx,Sy,W,M} <: BalanceLaw + slope_x::Sx + slope_y::Sy + width::W + mannings::M +end + +function OverlandFlowModel( + slope_x::Function, + slope_y::Function, + width::Function; + mannings::Function = (x, y) -> convert(eltype(x), 0.03)) + args = ( + slope_x, + slope_y, + width, + 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(area::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 + +function flux_first_order!(land::LandModel, surface::NoSurfaceFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +end + +""" + flux_first_order!(land::LandModel, surface::OverlandFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) + +The first order flux method for the OverlandFlow model. +""" +function flux_first_order!(land::LandModel, surface::OverlandFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) + x = aux.x + y = aux.y + width = surface.width(x, y) + area = max(eltype(state)(0.0), state.surface.area) + height = area / width + v = calculate_velocity(surface, x, y, height) + Q = area * v + flux.surface.area = Q +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. +""" +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.area = bc_function(aux⁻, t) +end + + +end diff --git a/src/Land/Model/land_bc.jl b/src/Land/Model/land_bc.jl index d1699a78dc6..7ea2e6c09dc 100644 --- a/src/Land/Model/land_bc.jl +++ b/src/Land/Model/land_bc.jl @@ -119,7 +119,7 @@ Base.@kwdef struct LandDomainBC{TBC, BBC, MinXBC, MaxXBC, MinYBC, MaxYBC} end """ - LandComponentBC{SW, SH} + LandComponentBC{SW, SH, SF} An object that holds the boundary conditions for each of the subcomponents of the land model. @@ -131,11 +131,11 @@ needs to define the BC for the components they wish to model. Base.@kwdef struct LandComponentBC{ SW <: AbstractBoundaryConditions, SH <: AbstractBoundaryConditions, - R <: AbstractBoundaryConditions, + SF <: AbstractBoundaryConditions, } soil_water::SW = NoBC() soil_heat::SH = NoBC() - river::R = NoBC() + surface::SF = NoBC() end @@ -184,7 +184,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...) - river_boundary_state!(nf, bc.river, land.river, land, args...) + surface_boundary_state!(nf, bc.surface, land.surface, land, args...) end @@ -222,5 +222,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...) - river_boundary_flux!(nf, bc.river, land.river, land, args...) + surface_boundary_flux!(nf, bc.surface, land.surface, land, args...) end diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index 5c80c079c1e..7c8de0645f4 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -62,7 +62,7 @@ function land_source!( t::Real, direction, ) - source.river.area += source_type(aux.x, aux.y, t) + source.surface.area += source_type(aux.x, aux.y, t) end function land_source!( diff --git a/test/Land/Model/Artifacts.toml b/test/Land/Model/Artifacts.toml index 1366f39123f..1dfef8738e8 100644 --- a/test/Land/Model/Artifacts.toml +++ b/test/Land/Model/Artifacts.toml @@ -1,2 +1,5 @@ [richards] git-tree-sha1 = "ff73fa6a0b6a807e71a6921f7ef7d0befe776edd" + +[tiltedv] +git-tree-sha1 = "db27235cb7ce2b7674607876da15d1635906b512" diff --git a/test/Land/Model/runtests.jl b/test/Land/Model/runtests.jl index e532b5321eb..a21199f2664 100644 --- a/test/Land/Model/runtests.jl +++ b/test/Land/Model/runtests.jl @@ -5,5 +5,5 @@ using Test, Pkg include("prescribed_twice.jl") include("freeze_thaw_alone.jl") include("test_runoff_functions.jl") - include("test_no_river.jl") + include("test_overland_flow.jl") end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 87634058a04..055074fe870 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -1,6 +1,3 @@ -# Test that the land model still runs, even with the lowest/simplest -# version of soil (prescribed heat and prescribed water - no state -# variables) using MPI using OrderedCollections using StaticArrays @@ -13,7 +10,7 @@ const param_set = EarthParameterSet() using ClimateMachine using ClimateMachine.Land -using ClimateMachine.Land.River +using ClimateMachine.Land.SurfaceFlow using ClimateMachine.Land.SoilWaterParameterizations using ClimateMachine.Mesh.Topologies using ClimateMachine.Mesh.Grids @@ -27,12 +24,11 @@ using ClimateMachine.VariableTemplates using ClimateMachine.SingleStackUtils using ClimateMachine.BalanceLaws: BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state +using ClimateMachine.ArtifactWrappers -function init_land_model!(land, state, aux, localgeo, time) -end - - -@testset "NoRiver Model" begin +@testset "NoSurfaceFlow Model" begin + function init_land_model!(land, state, aux, localgeo, time) + end ClimateMachine.init() FT = Float64 @@ -41,13 +37,13 @@ end soil_param_functions = nothing m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) - m_river = NoRiverModel() + m_surface = NoSurfaceFlowModel() sources = () m = LandModel( param_set, m_soil, - m_river; + m_surface; source = sources, init_state_prognostic = init_land_model!, ) @@ -97,23 +93,27 @@ end 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] == 3 @test size(Base.collect(keys(state_vars)))[1] == 0 end -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 +@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 -@testset "Analytical River 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 @@ -123,23 +123,20 @@ end m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model) - m_river = RiverModel( + m_surface = OverlandFlowModel( (x,y) -> eltype(x)(-0.0016), (x,y) -> eltype(x)(0.0), -# (x,y) -> eltype(x)(0.0016), (x,y) -> eltype(x)(1); mannings = (x,y) -> eltype(x)(0.025) ) - # Analytical test case defined as Model 1 in DOI: - # 10.1061/(ASCE)0733-9429(2007)133:2(217) - # Eqn 6 + bc = LandDomainBC( - minx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), + minx_bc = LandComponentBC(surface = Dirichlet((aux, t) -> eltype(aux)(0))), ) function init_land_model!(land, state, aux, localgeo, time) - state.river.area = eltype(state)(0) + state.surface.area = eltype(state)(0) end # units in m / s @@ -150,7 +147,7 @@ end m = LandModel( param_set, m_soil, - m_river; + m_surface; boundary_conditions = bc, source = sources, init_state_prognostic = init_land_model!, @@ -195,7 +192,7 @@ end Q = solver_config.Q area_index = - varsindex(vars_state(m, Prognostic(), FT), :river, :area) + varsindex(vars_state(m, Prognostic(), FT), :surface, :area) n_outputs = 60 every_x_simulation_time = ceil(Int, timeend / n_outputs) @@ -276,29 +273,39 @@ end end # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a GPU array) - @test sqrt_rmse_over_max_q = sqrt(mean((analytic.(time_data, alpha, t_c, t_r, i, L, m) .- q).^4.0))/ maximum(q) < 3e-3 + @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 -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)) +@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 - zout = zbase+zleft+zright - x, y, z = xin, yin, zout - return x, y, z -end - -#@testset "V Catchment Maxwell River Model" begin ClimateMachine.init() FT = Float64 @@ -329,7 +336,7 @@ end return x >= MFT(800) && x <= MFT(820) ? MFT(2.5*60 * 10^-3) : MFT(2.5*60 * 10^-4) end - m_river = RiverModel( + m_surface = OverlandFlowModel( x_slope, y_slope, (x,y) -> eltype(x)(1); @@ -337,13 +344,13 @@ end ) bc = LandDomainBC( - miny_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), - minx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), - maxx_bc = LandComponentBC(river = Dirichlet((aux, t) -> eltype(aux)(0))), + 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.river.area = eltype(state)(0) + state.surface.area = eltype(state)(0) end # units in m / s @@ -354,7 +361,7 @@ end m = LandModel( param_set, m_soil, - m_river; + m_surface; boundary_conditions = bc, source = sources, init_state_prognostic = init_land_model!, @@ -400,7 +407,7 @@ end Q = solver_config.Q area_index = - varsindex(vars_state(m, Prognostic(), FT), :river, :area) + varsindex(vars_state(m, Prognostic(), FT), :surface, :area) n_outputs = 60 every_x_simulation_time = ceil(Int, timeend / n_outputs) @@ -424,30 +431,26 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) - # Compare flowrate analytical derivation aux = solver_config.dg.state_auxiliary; x = aux[:,1,:] y = aux[:,2,:] z = aux[:,3,:] - # get all nodal points at the max X bound of the domain - mask2 = (Float64.(aux[:,2,:] .== 1000.0)) .==1# .+ Float64.(aux[:,1,:] .<= 820) .+ Float64.(aux[:,1,:] .> 800)) .== 3 - n_outputs = length(dons) + # 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# width = 1 here, everywhere - v = calculate_velocity(m_river,xv, 1000.0, height)# put in y = 1000.0 + 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 - # Maybe there is a conceptual error here? 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 - - - # get similation timesteps (s) - time_data = [dons[l]["t"][1] for l in 1:n_outputs] - #also helpful - #scatter(aux[:,1,:][:], aux[:,2,:][:], dons[50]["area"][:]) -#end +end diff --git a/test/Land/Model/test_overland_flow.jl b/test/Land/Model/test_overland_flow.jl new file mode 100644 index 00000000000..055074fe870 --- /dev/null +++ b/test/Land/Model/test_overland_flow.jl @@ -0,0 +1,456 @@ +using MPI +using OrderedCollections +using StaticArrays +using Test +using Statistics + +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 + +@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, + 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(60) + 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 + + +@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, 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), + (x,y) -> eltype(x)(1); + 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.area = 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, + 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 + + 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,)) + + # 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 area from nodal state (m^2) + area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] + height = area ./ ymax + # 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 + + # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a 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, + 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 From afd5128994e346a67c074e447563d858b5a6bfb2 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 11:31:40 -0700 Subject: [PATCH 15/24] soil tests seem to pass --- docs/src/APIs/Land/LandModel.md | 1 + src/Land/Model/LandModel.jl | 14 ++++++-------- src/Land/Model/prog_types.jl | 1 + src/Land/Model/prognostic_vars.jl | 7 ++++++- src/Land/Model/source.jl | 28 ++++++++++++++++----------- test/Land/Model/test_overland_flow.jl | 1 + 6 files changed, 32 insertions(+), 20 deletions(-) diff --git a/docs/src/APIs/Land/LandModel.md b/docs/src/APIs/Land/LandModel.md index 0f98b01fd6e..555e1bf4fb9 100644 --- a/docs/src/APIs/Land/LandModel.md +++ b/docs/src/APIs/Land/LandModel.md @@ -22,6 +22,7 @@ ClimateMachine.Land.SoilParamFunctions ClimateMachine.Land.get_water_content ClimateMachine.Land.get_temperature ClimateMachine.Land.PhaseChange +ClimateMachine.Land.Precip ``` ## Boundary Conditions diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 5b7c0631aa5..5115c1360ac 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -77,8 +77,8 @@ parameter_set(m::LandModel) = m.param_set """ LandModel( param_set::AbstractParameterSet, - soil::BalanceLaw, - surface::BalanceLaw; + soil::BalanceLaw; + surface::BalanceLaw = NoSurfaceFlowModel(), boundary_conditions::LBC = (), source::SRC = (), init_state_prognostic::IS = nothing @@ -88,8 +88,8 @@ Constructor for the LandModel structure. """ function LandModel( param_set::AbstractParameterSet, - soil::BalanceLaw, - surface::BalanceLaw; + soil::BalanceLaw; + surface::BalanceLaw = NoSurfaceFlowModel(), boundary_conditions::LBC = LandDomainBC(), source::SRC = (), init_state_prognostic::IS = nothing, @@ -175,7 +175,6 @@ function compute_gradient_argument!( ) compute_gradient_argument!(land, land.soil, transform, state, aux, t) - compute_gradient_argument!(land, land.surface, transform, state, aux, t) end function compute_gradient_flux!( @@ -279,13 +278,12 @@ 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("land_tendencies.jl") include("source.jl") -include("SurfaceFlow.jl") -using .SurfaceFlow function wavespeed( m::LandModel, diff --git a/src/Land/Model/prog_types.jl b/src/Land/Model/prog_types.jl index 77cab162a04..e1af7cc2f3b 100644 --- a/src/Land/Model/prog_types.jl +++ b/src/Land/Model/prog_types.jl @@ -5,3 +5,4 @@ struct VolumetricLiquidFraction <: AbstractPrognosticVariable end struct VolumetricInternalEnergy <: AbstractPrognosticVariable end struct VolumetricIceFraction <: AbstractPrognosticVariable end +struct SurfaceWaterArea <: AbstractPrognosticVariable end diff --git a/src/Land/Model/prognostic_vars.jl b/src/Land/Model/prognostic_vars.jl index acd841e4cf8..5d039b54b64 100644 --- a/src/Land/Model/prognostic_vars.jl +++ b/src/Land/Model/prognostic_vars.jl @@ -7,10 +7,15 @@ prognostic_vars(water::SoilWaterModel) = (VolumetricLiquidFraction(), VolumetricIceFraction()) prognostic_vars(heat::PrescribedTemperatureModel) = () prognostic_vars(heat::SoilHeatModel) = (VolumetricInternalEnergy(),) +prognostic_vars(surface::NoSurfaceFlowModel) = () +prognostic_vars(surface::OverlandFlowModel) = (SurfaceWaterArea(),) prognostic_vars(land::LandModel) = - (prognostic_vars(land.soil.water)..., prognostic_vars(land.soil.heat)...) + (prognostic_vars(land.soil.water)..., prognostic_vars(land.soil.heat)...,prognostic_vars(land.surface)...) + + 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, ::SurfaceWaterArea) = (state.surface, :area) diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index 41c1d6e1e5c..022d27dae0e 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -31,8 +31,13 @@ function precompute(land::LandModel, args, tt::Source) return (; dtup) end -#TODO: precip docs -struct Precip{FT, F} <: LandSource{FT} +""" + 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} @@ -44,17 +49,18 @@ function (p::Precip{FT})(x, y, t) where {FT} FT(p.precip(x, y, t)) end -function land_source!( - source_type::Precip, + +prognostic_vars(::Precip) = + (SurfaceWaterArea()) + +function source( + ::SurfaceWaterArea, + s::Precip, land::LandModel, - source::Vars, - state::Vars, - diffusive::Vars, - aux::Vars, - t::Real, - direction, + args, ) - source.surface.area += source_type(aux.x, aux.y, t) + @unpack aux, t = args + return s(aux.x, aux.y, t) end function precompute(source_type::PhaseChange, land::LandModel, args, tt::Source) diff --git a/test/Land/Model/test_overland_flow.jl b/test/Land/Model/test_overland_flow.jl index 055074fe870..bd0b505017a 100644 --- a/test/Land/Model/test_overland_flow.jl +++ b/test/Land/Model/test_overland_flow.jl @@ -3,6 +3,7 @@ using OrderedCollections using StaticArrays using Test using Statistics +using DelimitedFiles using CLIMAParameters struct EarthParameterSet <: AbstractEarthParameterSet end From f4bb38a43ea5d630139768e88241afe3fd2b703d Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 12:40:11 -0700 Subject: [PATCH 16/24] overland seems to work --- docs/list_of_apis.jl | 1 + docs/src/APIs/Land/LandModel.md | 1 - docs/src/APIs/Land/SurfaceFlow.md | 3 + src/Land/Model/LandModel.jl | 26 +-- src/Land/Model/River.jl | 89 +++++---- src/Land/Model/SurfaceFlow.jl | 150 ++++++++++---- src/Land/Model/land_bc.jl | 9 +- src/Land/Model/land_tendencies.jl | 25 ++- src/Land/Model/prog_types.jl | 1 - src/Land/Model/prognostic_vars.jl | 7 +- src/Land/Model/source.jl | 6 +- test/Land/Model/test_no_river.jl | 258 ++++++++++++------------ test/Land/Model/test_overland_flow.jl | 270 ++++++++++++++------------ 13 files changed, 506 insertions(+), 340 deletions(-) diff --git a/docs/list_of_apis.jl b/docs/list_of_apis.jl index 2f56a536667..602e78d3c2d 100644 --- a/docs/list_of_apis.jl +++ b/docs/list_of_apis.jl @@ -23,6 +23,7 @@ apis = [ "Soil Heat Parameterizations" => "APIs/Land/SoilHeatParameterizations.md", "Surface Flow" => "APIs/Land/SurfaceFlow.md", + "Radiative Boundary Conditions" => "RadiativeEnergyFlux.md", ], "Common" => [ "Orientations" => "APIs/Common/Orientations.md", diff --git a/docs/src/APIs/Land/LandModel.md b/docs/src/APIs/Land/LandModel.md index 555e1bf4fb9..0f98b01fd6e 100644 --- a/docs/src/APIs/Land/LandModel.md +++ b/docs/src/APIs/Land/LandModel.md @@ -22,7 +22,6 @@ ClimateMachine.Land.SoilParamFunctions ClimateMachine.Land.get_water_content ClimateMachine.Land.get_temperature ClimateMachine.Land.PhaseChange -ClimateMachine.Land.Precip ``` ## Boundary Conditions diff --git a/docs/src/APIs/Land/SurfaceFlow.md b/docs/src/APIs/Land/SurfaceFlow.md index ac8beca7920..78e470663e3 100644 --- a/docs/src/APIs/Land/SurfaceFlow.md +++ b/docs/src/APIs/Land/SurfaceFlow.md @@ -11,4 +11,7 @@ NoSurfaceFlowModel surface_boundary_flux! surface_boundary_state! calculate_velocity +Precip +VolumeAdvection +SurfaceWaterArea ``` \ No newline at end of file diff --git a/src/Land/Model/LandModel.jl b/src/Land/Model/LandModel.jl index 5115c1360ac..e335dedab2d 100644 --- a/src/Land/Model/LandModel.jl +++ b/src/Land/Model/LandModel.jl @@ -160,9 +160,18 @@ function flux_first_order!( state::Vars, aux::Vars, t::Real, - directions, -) - flux_first_order!(land, land.surface, flux, state, aux, t, directions) + 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 @@ -285,17 +294,10 @@ include("prognostic_vars.jl") include("land_tendencies.jl") include("source.jl") -function wavespeed( - m::LandModel, - n⁻, - state::Vars, - aux::Vars, - t::Real, - direction -) +function wavespeed(m::LandModel, n⁻, state::Vars, aux::Vars, t::Real, direction) FT = eltype(state) g = FT(9.8) - width = m.surface.width(aux.x,aux.y) + width = m.surface.width(aux.x, aux.y) area = max(eltype(state)(0.0), state.surface.area) height = area ./ width v = calculate_velocity(m.surface, aux.x, aux.y, height) diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl index c7058c722eb..c53ac32add7 100644 --- a/src/Land/Model/River.jl +++ b/src/Land/Model/River.jl @@ -1,24 +1,28 @@ module River - + using DocStringExtensions using ..Land using ..VariableTemplates import ..BalanceLaws: BalanceLaw, vars_state, - flux_first_order!, + flux_first_order!, Prognostic, - Auxiliary, + Auxiliary, Gradient, GradientFlux using ...DGMethods: LocalGeometry using StaticArrays: SVector -export RiverModel, NoRiverModel, river_boundary_flux!, river_boundary_state!, calculate_velocity +export RiverModel, + NoRiverModel, + river_boundary_flux!, + river_boundary_state!, + calculate_velocity struct NoRiverModel <: BalanceLaw end -struct RiverModel{Sx,Sy,W,M} <: BalanceLaw +struct RiverModel{Sx, Sy, W, M} <: BalanceLaw slope_x::Sx slope_y::Sy width::W @@ -29,13 +33,9 @@ function RiverModel( slope_x::Function, slope_y::Function, width::Function; - mannings::Function = (x, y) -> convert(eltype(x), 0.03)) - args = ( - slope_x, - slope_y, - width, - mannings, - ) + mannings::Function = (x, y) -> convert(eltype(x), 0.03), +) + args = (slope_x, slope_y, width, mannings) return RiverModel{typeof.(args)...}(args...) end @@ -44,25 +44,51 @@ function calculate_velocity(river, x::Real, y::Real, h::Real) sx = FT(river.slope_x(x, y)) sy = FT(river.slope_y(x, y)) mannings_coeff = FT(river.mannings(x, y)) - coeff = h^FT(2/3) / mannings_coeff + 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)) + return SVector( + -sign(sx) * coeff * sqrt(abs(sx)), + -sign(sy) * coeff * sqrt(abs(sy)), + zero(FT), + ) end vars_state(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) -function Land.land_init_aux!(land::LandModel, river::Union{NoRiverModel,RiverModel}, aux, geom::LocalGeometry) -end - -function Land.land_nodal_update_auxiliary_state!(land::LandModel, river::Union{NoRiverModel,RiverModel}, state, aux, t) -end +function Land.land_init_aux!( + land::LandModel, + river::Union{NoRiverModel, RiverModel}, + aux, + geom::LocalGeometry, +) end -function flux_first_order!(land::LandModel, river::NoRiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) -end +function Land.land_nodal_update_auxiliary_state!( + land::LandModel, + river::Union{NoRiverModel, RiverModel}, + state, + aux, + t, +) end -function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +function flux_first_order!( + land::LandModel, + river::NoRiverModel, + flux::Grad, + state::Vars, + aux::Vars, + t::Real, + directions, +) end + +function flux_first_order!( + land::LandModel, + river::RiverModel, + flux::Grad, + state::Vars, + aux::Vars, + t::Real, + directions, +) x = aux.x y = aux.y width = river.width(x, y) @@ -71,7 +97,7 @@ function flux_first_order!(land::LandModel, river::RiverModel, flux::Grad, state v = calculate_velocity(river, x, y, height) Q = area * v flux.river.area = Q -end +end # Boundary Conditions @@ -79,20 +105,18 @@ end function river_boundary_flux!( nf, bc::Land.AbstractBoundaryConditions, - m::Union{NoRiverModel,RiverModel}, + m::Union{NoRiverModel, RiverModel}, land::LandModel, _..., -) -end +) end function river_boundary_state!( nf, bc::Land.AbstractBoundaryConditions, - m::Union{NoRiverModel,RiverModel}, + m::Union{NoRiverModel, RiverModel}, land::LandModel, _..., -) -end +) end # Dirichlet BC for River function river_boundary_flux!( @@ -107,8 +131,7 @@ function river_boundary_flux!( aux⁻::Vars, t, _..., -) -end +) end function river_boundary_state!( nf, diff --git a/src/Land/Model/SurfaceFlow.jl b/src/Land/Model/SurfaceFlow.jl index 1b9f6a119bf..1ce41509ab8 100644 --- a/src/Land/Model/SurfaceFlow.jl +++ b/src/Land/Model/SurfaceFlow.jl @@ -1,21 +1,43 @@ module SurfaceFlow - + using DocStringExtensions +using UnPack using ..Land using ..VariableTemplates +using ..BalanceLaws import ..BalanceLaws: BalanceLaw, + prognostic_vars, + flux, + source, + precompute, + eq_tends, vars_state, - flux_first_order!, Prognostic, - Auxiliary, + Auxiliary, Gradient, GradientFlux -using ...DGMethods: LocalGeometry +using ...DGMethods: LocalGeometry, DGModel using StaticArrays: SVector -export OverlandFlowModel, NoSurfaceFlowModel, surface_boundary_flux!, surface_boundary_state!, calculate_velocity +export OverlandFlowModel, + NoSurfaceFlowModel, + surface_boundary_flux!, + surface_boundary_state!, + calculate_velocity, + Precip, + VolumeAdvection, + SurfaceWaterArea + +""" + SurfaceWaterArea <: AbstractPrognosticVariable + +The prognostic variable type for the overland flow model. Used only for +dispatching on. +""" +struct SurfaceWaterArea <: AbstractPrognosticVariable end + """ NoSurfaceFlowModel <: BalanceLaw @@ -33,7 +55,7 @@ to the LandModel reflecting the surface area of water per unit length. This model simulates depth averaged shallow water equation under the kinematic approximation, along with Manning's relationship. """ -struct OverlandFlowModel{Sx,Sy,W,M} <: BalanceLaw +struct OverlandFlowModel{Sx, Sy, W, M} <: BalanceLaw slope_x::Sx slope_y::Sy width::W @@ -44,13 +66,9 @@ function OverlandFlowModel( slope_x::Function, slope_y::Function, width::Function; - mannings::Function = (x, y) -> convert(eltype(x), 0.03)) - args = ( - slope_x, - slope_y, - width, - mannings, - ) + mannings::Function = (x, y) -> convert(eltype(x), 0.03), +) + args = (slope_x, slope_y, width, mannings) return OverlandFlowModel{typeof.(args)...}(args...) end @@ -66,39 +84,68 @@ function calculate_velocity(surface, x::Real, y::Real, h::Real) 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 + 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)) + 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(area::FT) -function Land.land_init_aux!(land::LandModel, surface::Union{NoSurfaceFlowModel,OverlandFlowModel}, aux, geom::LocalGeometry) -end +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 +function Land.land_nodal_update_auxiliary_state!( + land::LandModel, + surface::Union{NoSurfaceFlowModel, OverlandFlowModel}, + state, + aux, + t, +) end -function flux_first_order!(land::LandModel, surface::NoSurfaceFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) -end +""" + VolumeAdvection <: TendencyDef{Flux{FirstOrder}} +A first order flux type for the overland flow model. """ - flux_first_order!(land::LandModel, surface::OverlandFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +struct VolumeAdvection <: TendencyDef{Flux{FirstOrder}} end + +#function flux_first_order!(land::LandModel, surface::NoSurfaceFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) +#end -The first order flux method for the OverlandFlow model. """ -function flux_first_order!(land::LandModel, surface::OverlandFlowModel, flux::Grad, state::Vars, aux::Vars, t::Real, directions) + flux(::SurfaceWaterArea, ::VolumeAdvection, land::LandModel, args,) + +A first order flux method for the OverlandFlow model, adding in advection of water volume due to flowing water. +""" +function flux(::SurfaceWaterArea, ::VolumeAdvection, land::LandModel, args) + @unpack state, aux = args x = aux.x y = aux.y - width = surface.width(x, y) + width = land.surface.width(x, y) area = max(eltype(state)(0.0), state.surface.area) height = area / width - v = calculate_velocity(surface, x, y, height) + v = calculate_velocity(land.surface, x, y, height) Q = area * v - flux.surface.area = Q -end + return Q +end +#function flux_first_order!(land::LandModel, surface::OverlandFlowModel, fl#ux::Grad, state::Vars, aux::Vars, t::Real, directions) +# x = aux.x +# y = aux.y +# width = surface.width(x, y) +# area = max(eltype(state)(0.0), state.surface.area) +# height = area / width +# v = calculate_velocity(surface, x, y, height) +# Q = area * v +# flux.surface.area = Q +#end # Boundary Conditions @@ -106,20 +153,18 @@ end function surface_boundary_flux!( nf, bc::Land.AbstractBoundaryConditions, - m::Union{NoSurfaceFlowModel,OverlandFlowModel}, + m::Union{NoSurfaceFlowModel, OverlandFlowModel}, land::LandModel, _..., -) -end +) end function surface_boundary_state!( nf, bc::Land.AbstractBoundaryConditions, - m::Union{NoSurfaceFlowModel,OverlandFlowModel}, + m::Union{NoSurfaceFlowModel, OverlandFlowModel}, land::LandModel, _..., -) -end +) end """ surface_boundary_flux!( @@ -140,8 +185,7 @@ function surface_boundary_flux!( model::OverlandFlowModel, land::LandModel, _..., -) -end +) end """ surface_boundary_state!( @@ -178,5 +222,35 @@ function surface_boundary_state!( state⁺.surface.area = 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) = (SurfaceWaterArea(),) +#I thought there was a default so I wouldnt have to define this? +function precompute(source_type::Precip, land::LandModel, args, tt::Source) + nothing +end + + +function source(::SurfaceWaterArea, 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 dabc33dfb87..530f5b53e91 100644 --- a/src/Land/Model/land_bc.jl +++ b/src/Land/Model/land_bc.jl @@ -182,7 +182,14 @@ faces, as defined in the Driver configuration. """ function boundary_conditions(land::LandModel) bc = land.boundary_conditions - mytuple = (bc.bottom_bc, bc.surface_bc, bc.minx_bc, bc.maxx_bc, bc.miny_bc, bc.maxy_bc) + 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 diff --git a/src/Land/Model/land_tendencies.jl b/src/Land/Model/land_tendencies.jl index cd4b0892dc8..e7e7934f245 100644 --- a/src/Land/Model/land_tendencies.jl +++ b/src/Land/Model/land_tendencies.jl @@ -9,19 +9,36 @@ 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)..., +) +# can we replace these empty ones with a method for ::BalanceLaw? +eq_tends(::PV, ::AbstractSoilComponentModel, ::Flux{FirstOrder}) where {PV} = () +eq_tends(::PV, ::NoSurfaceFlowModel, ::Flux{FirstOrder}) where {PV} = () + +eq_tends(::SurfaceWaterArea, ::OverlandFlowModel, ::Flux{FirstOrder}) = + (VolumeAdvection(),) ##### ##### Second order fluxes ##### # Empty by default + +# can we replace these empty ones with a method for ::BalanceLaw? 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 +47,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 +54,6 @@ eq_tends( ::Flux{SecondOrder}, ) = (DiffHeatFlux(), DarcyDrivenHeatFlux()) -# TODO: move to soi_water.jl? eq_tends( ::VolumetricInternalEnergy, ::SoilHeatModel, diff --git a/src/Land/Model/prog_types.jl b/src/Land/Model/prog_types.jl index e1af7cc2f3b..77cab162a04 100644 --- a/src/Land/Model/prog_types.jl +++ b/src/Land/Model/prog_types.jl @@ -5,4 +5,3 @@ struct VolumetricLiquidFraction <: AbstractPrognosticVariable end struct VolumetricInternalEnergy <: AbstractPrognosticVariable end struct VolumetricIceFraction <: AbstractPrognosticVariable end -struct SurfaceWaterArea <: AbstractPrognosticVariable end diff --git a/src/Land/Model/prognostic_vars.jl b/src/Land/Model/prognostic_vars.jl index 5d039b54b64..e9888852090 100644 --- a/src/Land/Model/prognostic_vars.jl +++ b/src/Land/Model/prognostic_vars.jl @@ -10,8 +10,11 @@ prognostic_vars(heat::SoilHeatModel) = (VolumetricInternalEnergy(),) prognostic_vars(surface::NoSurfaceFlowModel) = () prognostic_vars(surface::OverlandFlowModel) = (SurfaceWaterArea(),) -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)..., + prognostic_vars(land.surface)..., +) diff --git a/src/Land/Model/source.jl b/src/Land/Model/source.jl index 022d27dae0e..f9631888567 100644 --- a/src/Land/Model/source.jl +++ b/src/Land/Model/source.jl @@ -1,5 +1,5 @@ #### Land sources -export PhaseChange, Precip +export PhaseChange#, Precip function heaviside(x::FT) where {FT} if x >= FT(0) @@ -30,7 +30,7 @@ function precompute(land::LandModel, args, tt::Source) end) return (; dtup) end - +#= """ Precip <: TendencyDef{Source} @@ -62,7 +62,7 @@ function source( @unpack aux, t = args return s(aux.x, aux.y, t) end - +=# function precompute(source_type::PhaseChange, land::LandModel, args, tt::Source) @unpack state, diffusive, aux, t, direction = args diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl index 055074fe870..d9484de3d65 100644 --- a/test/Land/Model/test_no_river.jl +++ b/test/Land/Model/test_no_river.jl @@ -27,8 +27,7 @@ using ClimateMachine.BalanceLaws: using ClimateMachine.ArtifactWrappers @testset "NoSurfaceFlow Model" begin - function init_land_model!(land, state, aux, localgeo, time) - end + function init_land_model!(land, state, aux, localgeo, time) end ClimateMachine.init() FT = Float64 @@ -105,10 +104,17 @@ end # 10.1061/(ASCE)0733-9429(2007)133:2(217) # Eqn 6 - - function warp_constant_slope(xin, yin, zin; topo_max = 0.2, zmin = -0.1, xmax = 400) + + 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) + 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 @@ -122,19 +128,21 @@ end 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), - (x,y) -> eltype(x)(1); - mannings = (x,y) -> eltype(x)(0.025) + (x, y) -> eltype(x)(-0.0016), + (x, y) -> eltype(x)(0.0), + (x, y) -> eltype(x)(1); + mannings = (x, y) -> eltype(x)(0.025), ) - + bc = LandDomainBC( - minx_bc = LandComponentBC(surface = Dirichlet((aux, t) -> eltype(aux)(0))), + minx_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 @@ -153,33 +161,37 @@ end init_state_prognostic = init_land_model!, ) - N_poly = 1; + 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); + zmax = FT(0) + zmin = FT(-0.1) xmax = FT(182.88) ymax = FT(1.0) - topo_max = FT(0.0016*xmax) + topo_max = FT(0.0016 * xmax) driver_config = ClimateMachine.MultiColumnLandModel( "LandModel", (N_poly, N_poly), - (xres,yres,zres), + (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), - ); + meshwarp = (x...) -> warp_constant_slope( + x...; + topo_max = topo_max, + zmin = zmin, + xmax = xmax, + ), + ) t0 = FT(0) - timeend = FT(60*60) + timeend = FT(60 * 60) dt = FT(10) solver_config = ClimateMachine.SolverConfiguration( @@ -190,9 +202,8 @@ end ) mygrid = solver_config.dg.grid Q = solver_config.Q - - area_index = - varsindex(vars_state(m, Prognostic(), FT), :surface, :area) + + area_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :area) n_outputs = 60 every_x_simulation_time = ceil(Int, timeend / n_outputs) @@ -205,10 +216,7 @@ end ) do (init = false) t = ODESolvers.gettime(solver_config.solver) area = Q[:, area_index, :] - all_vars = Dict{String, Array}( - "t" => [t], - "area" => area, - ) + all_vars = Dict{String, Array}("t" => [t], "area" => area) dons[iostep[1]] = all_vars iostep[1] += 1 return @@ -217,63 +225,68 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) # Compare flowrate analytical derivation - aux = solver_config.dg.state_auxiliary; - + aux = solver_config.dg.state_auxiliary + # get all nodal points at the max X bound of the domain - mask = Array(aux[:,1,:] .== 182.88) + mask = Array(aux[:, 1, :] .== 182.88) n_outputs = length(dons) # get prognostic variable area from nodal state (m^2) area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] height = area ./ ymax # get similation timesteps (s) time_data = [dons[l]["t"][1] for l in 1:n_outputs] - - alpha = sqrt(0.0016)/0.025 + + 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) + 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) + + 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) + + function analytic(t, alpha, t_c, t_r, i, L, m) if t < t_c - return alpha*(i*t)^(m) + return alpha * (i * t)^(m) end - + if t <= t_r && t > t_c - return alpha*(i*t_c)^(m) + return alpha * (i * t_c)^(m) end - + if t > t_r - yL = (i*(t-t_r)) + yL = (i * (t - t_r)) delta = 1 - error = g(m,yL,i,t_r,L,alpha,t) + 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) + 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 - + return alpha * yL^m + end - + end - + # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a 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 + @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 @@ -288,76 +301,83 @@ end ),], ) 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)) + zbase = slope_v * (yin - FT(1000)) zleft = FT(0.0) zright = FT(0.0) if xin < FT(800) - zleft = slope_sides*(xin-FT(800)) + zleft = slope_sides * (xin - FT(800)) end if xin > FT(820) - zright = slope_sides*(FT(1)+(xin-FT(820))/FT(800)) + zright = slope_sides * (FT(1) + (xin - FT(820)) / FT(800)) end - zout = zbase+zleft+zright + 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 + 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) + 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); + (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))), + 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, @@ -366,37 +386,37 @@ end source = sources, init_state_prognostic = init_land_model!, ) - - N_poly_hori = 1; - N_poly_vert = 1; + + 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); + 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), + (xres, yres, zres), xmax, ymax, zmax, param_set, m; zmin = zmin, - numerical_flux_first_order = RusanovNumericalFlux() + numerical_flux_first_order = RusanovNumericalFlux(), # meshwarp = (x...) -> warp_tilted_v(x...), - ); - + ) + t0 = FT(0) - timeend = FT(180*60) + timeend = FT(180 * 60) dt = FT(0.5) - + solver_config = ClimateMachine.SolverConfiguration( t0, timeend, @@ -405,49 +425,49 @@ end ) mygrid = solver_config.dg.grid Q = solver_config.Q - - area_index = - varsindex(vars_state(m, Prognostic(), FT), :surface, :area) + + 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, - ) + 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,:] + + 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 + 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 + 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 + 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)) diff --git a/test/Land/Model/test_overland_flow.jl b/test/Land/Model/test_overland_flow.jl index bd0b505017a..fb76164a8ae 100644 --- a/test/Land/Model/test_overland_flow.jl +++ b/test/Land/Model/test_overland_flow.jl @@ -28,8 +28,7 @@ using ClimateMachine.BalanceLaws: using ClimateMachine.ArtifactWrappers @testset "NoSurfaceFlow Model" begin - function init_land_model!(land, state, aux, localgeo, time) - end + function init_land_model!(land, state, aux, localgeo, time) end ClimateMachine.init() FT = Float64 @@ -43,8 +42,8 @@ using ClimateMachine.ArtifactWrappers sources = () m = LandModel( param_set, - m_soil, - m_surface; + m_soil; + surface = m_surface, source = sources, init_state_prognostic = init_land_model!, ) @@ -106,10 +105,17 @@ end # 10.1061/(ASCE)0733-9429(2007)133:2(217) # Eqn 6 - - function warp_constant_slope(xin, yin, zin; topo_max = 0.2, zmin = -0.1, xmax = 400) + + 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) + 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 @@ -123,19 +129,21 @@ end 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), - (x,y) -> eltype(x)(1); - mannings = (x,y) -> eltype(x)(0.025) + (x, y) -> eltype(x)(-0.0016), + (x, y) -> eltype(x)(0.0), + (x, y) -> eltype(x)(1); + mannings = (x, y) -> eltype(x)(0.025), ) - + bc = LandDomainBC( - minx_bc = LandComponentBC(surface = Dirichlet((aux, t) -> eltype(aux)(0))), + minx_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 @@ -147,40 +155,44 @@ end m = LandModel( param_set, - m_soil, - m_surface; + m_soil; + surface = m_surface, boundary_conditions = bc, source = sources, init_state_prognostic = init_land_model!, ) - N_poly = 1; + 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); + zmax = FT(0) + zmin = FT(-0.1) xmax = FT(182.88) ymax = FT(1.0) - topo_max = FT(0.0016*xmax) + topo_max = FT(0.0016 * xmax) driver_config = ClimateMachine.MultiColumnLandModel( "LandModel", (N_poly, N_poly), - (xres,yres,zres), + (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), - ); + meshwarp = (x...) -> warp_constant_slope( + x...; + topo_max = topo_max, + zmin = zmin, + xmax = xmax, + ), + ) t0 = FT(0) - timeend = FT(60*60) + timeend = FT(60 * 60) dt = FT(10) solver_config = ClimateMachine.SolverConfiguration( @@ -191,9 +203,8 @@ end ) mygrid = solver_config.dg.grid Q = solver_config.Q - - area_index = - varsindex(vars_state(m, Prognostic(), FT), :surface, :area) + + area_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :area) n_outputs = 60 every_x_simulation_time = ceil(Int, timeend / n_outputs) @@ -206,10 +217,7 @@ end ) do (init = false) t = ODESolvers.gettime(solver_config.solver) area = Q[:, area_index, :] - all_vars = Dict{String, Array}( - "t" => [t], - "area" => area, - ) + all_vars = Dict{String, Array}("t" => [t], "area" => area) dons[iostep[1]] = all_vars iostep[1] += 1 return @@ -218,63 +226,68 @@ end ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) # Compare flowrate analytical derivation - aux = solver_config.dg.state_auxiliary; - + aux = solver_config.dg.state_auxiliary + # get all nodal points at the max X bound of the domain - mask = Array(aux[:,1,:] .== 182.88) + mask = Array(aux[:, 1, :] .== 182.88) n_outputs = length(dons) # get prognostic variable area from nodal state (m^2) area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] height = area ./ ymax # get similation timesteps (s) time_data = [dons[l]["t"][1] for l in 1:n_outputs] - - alpha = sqrt(0.0016)/0.025 + + 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) + 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) + + 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) + + function analytic(t, alpha, t_c, t_r, i, L, m) if t < t_c - return alpha*(i*t)^(m) + return alpha * (i * t)^(m) end - + if t <= t_r && t > t_c - return alpha*(i*t_c)^(m) + return alpha * (i * t_c)^(m) end - + if t > t_r - yL = (i*(t-t_r)) + yL = (i * (t - t_r)) delta = 1 - error = g(m,yL,i,t_r,L,alpha,t) + 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) + 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 - + return alpha * yL^m + end - + end - + # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a 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 + @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 @@ -289,115 +302,122 @@ end ),], ) 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)) + zbase = slope_v * (yin - FT(1000)) zleft = FT(0.0) zright = FT(0.0) if xin < FT(800) - zleft = slope_sides*(xin-FT(800)) + zleft = slope_sides * (xin - FT(800)) end if xin > FT(820) - zright = slope_sides*(FT(1)+(xin-FT(820))/FT(800)) + zright = slope_sides * (FT(1) + (xin - FT(820)) / FT(800)) end - zout = zbase+zleft+zright + 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 + 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) + 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); + (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))), + 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, - m_surface; + m_soil; + surface = m_surface, boundary_conditions = bc, source = sources, init_state_prognostic = init_land_model!, ) - - N_poly_hori = 1; - N_poly_vert = 1; + + 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); + 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), + (xres, yres, zres), xmax, ymax, zmax, param_set, m; zmin = zmin, - numerical_flux_first_order = RusanovNumericalFlux() + numerical_flux_first_order = RusanovNumericalFlux(), # meshwarp = (x...) -> warp_tilted_v(x...), - ); - + ) + t0 = FT(0) - timeend = FT(180*60) + timeend = FT(180 * 60) dt = FT(0.5) - + solver_config = ClimateMachine.SolverConfiguration( t0, timeend, @@ -406,49 +426,49 @@ end ) mygrid = solver_config.dg.grid Q = solver_config.Q - - area_index = - varsindex(vars_state(m, Prognostic(), FT), :surface, :area) + + 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, - ) + 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,:] + + 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 + 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 + 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 + 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)) From 477f4e5a1a8d4ffeeef3550ddc94504e12cb35f5 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 12:50:38 -0700 Subject: [PATCH 17/24] deleted old River.jl and test_no_river.jl files from remote --- src/Land/Model/River.jl | 154 ---------- test/Land/Model/test_no_river.jl | 476 ------------------------------- 2 files changed, 630 deletions(-) delete mode 100644 src/Land/Model/River.jl delete mode 100644 test/Land/Model/test_no_river.jl diff --git a/src/Land/Model/River.jl b/src/Land/Model/River.jl deleted file mode 100644 index c53ac32add7..00000000000 --- a/src/Land/Model/River.jl +++ /dev/null @@ -1,154 +0,0 @@ -module River - -using DocStringExtensions -using ..Land -using ..VariableTemplates -import ..BalanceLaws: - BalanceLaw, - vars_state, - flux_first_order!, - Prognostic, - Auxiliary, - Gradient, - GradientFlux - -using ...DGMethods: LocalGeometry -using StaticArrays: SVector - -export RiverModel, - NoRiverModel, - river_boundary_flux!, - river_boundary_state!, - calculate_velocity - -struct NoRiverModel <: BalanceLaw end -struct RiverModel{Sx, Sy, W, M} <: BalanceLaw - slope_x::Sx - slope_y::Sy - width::W - mannings::M -end - -function RiverModel( - slope_x::Function, - slope_y::Function, - width::Function; - mannings::Function = (x, y) -> convert(eltype(x), 0.03), -) - args = (slope_x, slope_y, width, mannings) - return RiverModel{typeof.(args)...}(args...) -end - -function calculate_velocity(river, x::Real, y::Real, h::Real) - FT = eltype(h) - sx = FT(river.slope_x(x, y)) - sy = FT(river.slope_y(x, y)) - mannings_coeff = FT(river.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(river::RiverModel, st::Prognostic, FT) = @vars(area::FT) - -function Land.land_init_aux!( - land::LandModel, - river::Union{NoRiverModel, RiverModel}, - aux, - geom::LocalGeometry, -) end - -function Land.land_nodal_update_auxiliary_state!( - land::LandModel, - river::Union{NoRiverModel, RiverModel}, - state, - aux, - t, -) end - -function flux_first_order!( - land::LandModel, - river::NoRiverModel, - flux::Grad, - state::Vars, - aux::Vars, - t::Real, - directions, -) end - -function flux_first_order!( - land::LandModel, - river::RiverModel, - flux::Grad, - state::Vars, - aux::Vars, - t::Real, - directions, -) - x = aux.x - y = aux.y - width = river.width(x, y) - area = max(eltype(state)(0.0), state.river.area) - height = area / width - v = calculate_velocity(river, x, y, height) - Q = area * v - flux.river.area = Q -end - -# Boundary Conditions - -# General case - to be used with bc::NoBC -function river_boundary_flux!( - nf, - bc::Land.AbstractBoundaryConditions, - m::Union{NoRiverModel, RiverModel}, - land::LandModel, - _..., -) end - -function river_boundary_state!( - nf, - bc::Land.AbstractBoundaryConditions, - m::Union{NoRiverModel, RiverModel}, - land::LandModel, - _..., -) end - -# Dirichlet BC for River -function river_boundary_flux!( - nf, - bc::Land.Dirichlet, - model::RiverModel, - land::LandModel, - state⁺::Vars, - aux⁺::Vars, - nM, - state⁻::Vars, - aux⁻::Vars, - t, - _..., -) end - -function river_boundary_state!( - nf, - bc::Land.Dirichlet, - model::RiverModel, - land::LandModel, - state⁺::Vars, - aux⁺::Vars, - nM, - state⁻::Vars, - aux⁻::Vars, - t, - _..., -) - bc_function = bc.state_bc - state⁺.river.area = bc_function(aux⁻, t) -end - - -end diff --git a/test/Land/Model/test_no_river.jl b/test/Land/Model/test_no_river.jl deleted file mode 100644 index d9484de3d65..00000000000 --- a/test/Land/Model/test_no_river.jl +++ /dev/null @@ -1,476 +0,0 @@ -using MPI -using OrderedCollections -using StaticArrays -using Test -using Statistics - -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 - -@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, - 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(60) - 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 - - -@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, - 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), - (x, y) -> eltype(x)(1); - 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.area = 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, - 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 - - 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,)) - - # 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 area from nodal state (m^2) - area = [mean(Array(dons[k]["area"])[mask[:]]) for k in 1:n_outputs] - height = area ./ ymax - # 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 - - # The Ref's here are to ensure it works on CPU and GPU compatible array backends (q can be a 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, - 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 From f2f3dc40aa35f82b2788dbd8c946710e54354831 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 13:11:26 -0700 Subject: [PATCH 18/24] fixed unit test --- test/Land/Model/prescribed_twice.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Land/Model/prescribed_twice.jl b/test/Land/Model/prescribed_twice.jl index fd0ff045ff2..49bb4121fbd 100644 --- a/test/Land/Model/prescribed_twice.jl +++ b/test/Land/Model/prescribed_twice.jl @@ -93,6 +93,6 @@ using ClimateMachine.BalanceLaws: ) #Make sure it runs, and that there are no state variables, and only "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 From cb292993a25f9457116a57658bfa9551b87b20e9 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 13:34:11 -0700 Subject: [PATCH 19/24] removed old file, fixed doc --- docs/list_of_apis.jl | 2 +- test/Land/Model/runtests.jl | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/list_of_apis.jl b/docs/list_of_apis.jl index 602e78d3c2d..8125fc4ee0b 100644 --- a/docs/list_of_apis.jl +++ b/docs/list_of_apis.jl @@ -23,7 +23,7 @@ apis = [ "Soil Heat Parameterizations" => "APIs/Land/SoilHeatParameterizations.md", "Surface Flow" => "APIs/Land/SurfaceFlow.md", - "Radiative Boundary Conditions" => "RadiativeEnergyFlux.md", + "Radiative Boundary Conditions" => "APIs/Land/RadiativeEnergyFlux.md", ], "Common" => [ "Orientations" => "APIs/Common/Orientations.md", diff --git a/test/Land/Model/runtests.jl b/test/Land/Model/runtests.jl index 1bed83d99ea..803ba0f3ed2 100644 --- a/test/Land/Model/runtests.jl +++ b/test/Land/Model/runtests.jl @@ -4,7 +4,6 @@ using Test, Pkg include("test_water_parameterizations.jl") include("prescribed_twice.jl") include("freeze_thaw_alone.jl") - include("test_runoff_functions.jl") include("test_overland_flow.jl") include("test_physical_bc.jl") end From 422a50da5f85f8206781d889c3d3a0bc13775ee2 Mon Sep 17 00:00:00 2001 From: Katherine Deck Date: Thu, 8 Apr 2021 13:37:45 -0700 Subject: [PATCH 20/24] formatter doesnt do the right thing --- docs/list_of_apis.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/list_of_apis.jl b/docs/list_of_apis.jl index 8125fc4ee0b..d2fd7830cc8 100644 --- a/docs/list_of_apis.jl +++ b/docs/list_of_apis.jl @@ -23,7 +23,8 @@ apis = [ "Soil Heat Parameterizations" => "APIs/Land/SoilHeatParameterizations.md", "Surface Flow" => "APIs/Land/SurfaceFlow.md", - "Radiative Boundary Conditions" => "APIs/Land/RadiativeEnergyFlux.md", + "Radiative Boundary Conditions" => + "APIs/Land/RadiativeEnergyFlux.md", ], "Common" => [ "Orientations" => "APIs/Common/Orientations.md", From 5596ec898d76e432239777cbd1c9cd7c56315df6 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Fri, 9 Apr 2021 11:00:35 -0600 Subject: [PATCH 21/24] move long running tilted v into its own test, run on cluster --- .buildkite/pipeline.yml | 19 ++ test/Land/Model/Artifacts.toml | 5 +- test/Land/Model/runtests.jl | 2 +- ...flow.jl => test_overland_flow_analytic.jl} | 199 +---------------- .../Model/test_overland_flow_vcatchment.jl | 211 ++++++++++++++++++ 5 files changed, 239 insertions(+), 197 deletions(-) rename test/Land/Model/{test_overland_flow.jl => test_overland_flow_analytic.jl} (57%) create mode 100644 test/Land/Model/test_overland_flow_vcatchment.jl 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 From 8d015c5de189f414ab05c11f5ffcf763de767484 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Fri, 9 Apr 2021 12:57:04 -0600 Subject: [PATCH 22/24] fix boundry condition test to match new interface --- test/Land/Model/test_bc_3d.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/test/Land/Model/test_bc_3d.jl b/test/Land/Model/test_bc_3d.jl index 159da38c37b..9a0a45f3c20 100644 --- a/test/Land/Model/test_bc_3d.jl +++ b/test/Land/Model/test_bc_3d.jl @@ -44,13 +44,23 @@ 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) From a5ebc582dc72c1533644340b6c0c1cf4f7d264d6 Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Fri, 9 Apr 2021 13:59:14 -0600 Subject: [PATCH 23/24] gpu v catchment test fixes --- test/Land/Model/test_overland_flow_vcatchment.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/Land/Model/test_overland_flow_vcatchment.jl b/test/Land/Model/test_overland_flow_vcatchment.jl index c6b8a841571..d09b9a6ea27 100644 --- a/test/Land/Model/test_overland_flow_vcatchment.jl +++ b/test/Land/Model/test_overland_flow_vcatchment.jl @@ -183,11 +183,11 @@ using ClimateMachine.ArtifactWrappers ClimateMachine.invoke!(solver_config; user_callbacks = (callback,)) aux = solver_config.dg.state_auxiliary - x = aux[:, 1, :] - y = aux[:, 2, :] - z = aux[:, 3, :] + x = Array(aux[:, 1, :]) + y = Array(aux[:, 2, :]) + z = Array(aux[:, 3, :]) # Get points at outlet (y = ymax) - mask2 = (Float64.(aux[:, 2, :] .== 1000.0)) .== 1 + mask2 = (Float64.(y .== 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 @@ -207,5 +207,4 @@ using ClimateMachine.ArtifactWrappers ds_tv = readdlm(data, ',') error = sqrt(mean(Q .- ds_tv)) @test error < 1e-10 - end From 4053233c2b994d9c02f0d8a27ca82198710374fc Mon Sep 17 00:00:00 2001 From: Jake Bolewski Date: Mon, 12 Apr 2021 16:31:34 -0600 Subject: [PATCH 24/24] format --- test/Land/Model/test_bc_3d.jl | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/test/Land/Model/test_bc_3d.jl b/test/Land/Model/test_bc_3d.jl index 9a0a45f3c20..0464af9a85d 100644 --- a/test/Land/Model/test_bc_3d.jl +++ b/test/Land/Model/test_bc_3d.jl @@ -50,18 +50,10 @@ using ClimateMachine.BalanceLaws: bc = LandDomainBC( bottom_bc = LandComponentBC(soil_water = Neumann(bottom_flux)), surface_bc = LandComponentBC(soil_water = Dirichlet(surface_state)), - 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), - ), + 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()