From 89437d927844118bb6ef5c8d14a50853b4cf4ff3 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Fri, 13 Dec 2024 16:15:21 -0800 Subject: [PATCH] increased tf to 2 years --- .../standalone/Soil/evaporation_modified.jl | 34 ++++++++----------- experiments/long_runs/land.jl | 2 +- experiments/long_runs/soil.jl | 2 +- src/shared_utilities/implicit_timestepping.jl | 6 ++-- src/standalone/Soil/energy_hydrology.jl | 14 ++++---- 5 files changed, 27 insertions(+), 31 deletions(-) diff --git a/docs/tutorials/standalone/Soil/evaporation_modified.jl b/docs/tutorials/standalone/Soil/evaporation_modified.jl index a018b37f2a..eacea1be0e 100644 --- a/docs/tutorials/standalone/Soil/evaporation_modified.jl +++ b/docs/tutorials/standalone/Soil/evaporation_modified.jl @@ -62,7 +62,7 @@ e = rh * esat q = FT(0.622 * e / (101325 - 0.378 * e)) K_sat = FT(225.1 / 3600 / 24 / 1000) # precip = (t) -> min(-(K_sat/18) * sin(2*pi*t/(86400*3)), 0) -precip = (t) -> min((K_sat/0.1) * (sin(2*pi*t/(86400))+0.9), 0) +precip = (t) -> min((K_sat / 0.1) * (sin(2 * pi * t / (86400)) + 0.9), 0) T_atmos = (t) -> T_air u_atmos = (t) -> 1.0 q_atmos = (t) -> q @@ -246,19 +246,25 @@ data_e = ref_soln_E_350mm[:, 2]; outdir = joinpath(@__DIR__, "jacobian_van_g_evap_Ss1e-1") !ispath(outdir) && mkdir(outdir) S = [ - Array(parent(ClimaLand.top_center_to_surface((sol.u[k].soil.ϑ_l .- θ_r) ./ (ν .- θ_r))))[1] for k in 1:length(sol.t) + Array( + parent( + ClimaLand.top_center_to_surface( + (sol.u[k].soil.ϑ_l .- θ_r) ./ (ν .- θ_r), + ), + ), + )[1] for k in 1:length(sol.t) ]; θ_i = [ - Array(parent(ClimaLand.top_center_to_surface(sol.u[k].soil.θ_i)))[1] for k in 1:length(sol.t) + Array(parent(ClimaLand.top_center_to_surface(sol.u[k].soil.θ_i)))[1] for + k in 1:length(sol.t) ]; T = [ - Array(parent(ClimaLand.top_center_to_surface(sv.saveval[k].soil.T)))[1] for k in 1:length(sol.t) -]; -P_liq = [ - Array(parent(sv.saveval[k].drivers.P_liq))[1] for k in 1:length(sol.t) + Array(parent(ClimaLand.top_center_to_surface(sv.saveval[k].soil.T)))[1] + for k in 1:length(sol.t) ]; +P_liq = [Array(parent(sv.saveval[k].drivers.P_liq))[1] for k in 1:length(sol.t)]; fig = Figure(size = (800, 400)) ax = Axis( @@ -271,23 +277,13 @@ CairoMakie.lines!(ax, sol.t ./ 3600 ./ 24, FT.(S)) save(joinpath(outdir, "evaporation_lehmann2008_fig8b_water.png"), fig); fig = Figure(size = (800, 400)) -ax = Axis( - fig[1, 1], - xlabel = "Day", - ylabel = "temp", - title = "temp in soil", -) +ax = Axis(fig[1, 1], xlabel = "Day", ylabel = "temp", title = "temp in soil") CairoMakie.lines!(ax, sol.t ./ 3600 ./ 24, FT.(T)) save(joinpath(outdir, "evaporation_lehmann2008_fig8b_temp.png"), fig); fig = Figure(size = (800, 400)) -ax = Axis( - fig[1, 1], - xlabel = "Day", - ylabel = "precip/evap", - title = "precip", -) +ax = Axis(fig[1, 1], xlabel = "Day", ylabel = "precip/evap", title = "precip") CairoMakie.xlims!(minimum(data_dates), sol.t[end] ./ 3600 ./ 24) CairoMakie.lines!( ax, diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl index 505c97dbf0..e03f966150 100644 --- a/experiments/long_runs/land.jl +++ b/experiments/long_runs/land.jl @@ -379,7 +379,7 @@ end function setup_and_solve_problem(; greet = false) t0 = 0.0 - tf = 60 * 60.0 * 24 * 365 + tf = 60 * 60.0 * 24 * 365 * 2 Δt = 450.0 nelements = (101, 15) if greet diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl index e778799508..420f20bef8 100644 --- a/experiments/long_runs/soil.jl +++ b/experiments/long_runs/soil.jl @@ -220,7 +220,7 @@ end function setup_and_solve_problem(; greet = false) t0 = 0.0 - tf = 60 * 60.0 * 24 * 365 + tf = 60 * 60.0 * 24 * 365 * 2 Δt = 450.0 nelements = (101, 15) if greet diff --git a/src/shared_utilities/implicit_timestepping.jl b/src/shared_utilities/implicit_timestepping.jl index 59235c5b21..9a61edb933 100644 --- a/src/shared_utilities/implicit_timestepping.jl +++ b/src/shared_utilities/implicit_timestepping.jl @@ -153,8 +153,10 @@ function ImplicitEquationJacobian(Y::ClimaCore.Fields.FieldVector) # These are the off-diagonal terms in the Jacobian matrix # Here, we take the convention that each pair has order (T_x, y) to produce ∂T_x/∂y as above off_diagonal_pairs = ((@name(soil.ρe_int), @name(soil.ϑ_l)),) - available_off_diagonal_pairs = - MatrixFields.unrolled_filter(pair -> all(is_in_Y.(pair)),off_diagonal_pairs) + available_off_diagonal_pairs = MatrixFields.unrolled_filter( + pair -> all(is_in_Y.(pair)), + off_diagonal_pairs, + ) implicit_off_diagonals = MatrixFields.unrolled_map( pair -> (pair[1], pair[2]) => diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl index 6d7c482111..e7565387f5 100644 --- a/src/standalone/Soil/energy_hydrology.jl +++ b/src/standalone/Soil/energy_hydrology.jl @@ -322,7 +322,7 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} # The derivative of the residual with respect to the prognostic variable ∂ϑres∂ϑ = matrix[@name(soil.ϑ_l), @name(soil.ϑ_l)] ∂ρeres∂ρe = matrix[@name(soil.ρe_int), @name(soil.ρe_int)] - ∂ρeres∂ϑ = matrix[@name(soil.ρe_int), @name(soil.ϑ_l)] + ∂ρeres∂ϑ = matrix[@name(soil.ρe_int), @name(soil.ϑ_l)] # If the top BC is a `MoistureStateBC`, add the term from the top BC # flux before applying divergence if haskey(p.soil, :dfluxBCdY) @@ -376,16 +376,14 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} end @. ∂ρeres∂ϑ = -dtγ * ( - divf2c_matrix() ⋅ - MatrixFields.DiagonalMatrixRow( + divf2c_matrix() ⋅ MatrixFields.DiagonalMatrixRow( -interpc2f_op( volumetric_internal_energy_liq( p.soil.T, model.parameters.earth_param_set, ) * p.soil.K, - ) - ) ⋅ - gradc2f_matrix() ⋅ MatrixFields.DiagonalMatrixRow( + ), + ) ⋅ gradc2f_matrix() ⋅ MatrixFields.DiagonalMatrixRow( ClimaLand.Soil.dψdϑ( hydrology_cm, Y.soil.ϑ_l, @@ -395,9 +393,9 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} ), ) ) - (I,) - - + + @. ∂ρeres∂ρe = -dtγ * ( divf2c_matrix() ⋅