From a945721edd70aec1fdf91497ff08daf7664e9e3d Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 1 Oct 2024 20:56:02 +0200 Subject: [PATCH 1/2] Switch pow to FastPower --- Manifest.toml | 23 ++++++++++++++++++++++- Project.toml | 2 ++ src/Wflow.jl | 1 + src/utils.jl | 4 ---- test/horizontal_process.jl | 8 ++++---- test/reservoir_lake.jl | 16 ++++++++-------- test/run_sbm.jl | 2 +- test/vertical_process.jl | 4 ++-- 8 files changed, 40 insertions(+), 20 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index a72a3f92b..b4ebac40f 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "acf9296ffbaffd3f335bdb6005badef66b3e10fe" +project_hash = "6350b5664a0cf56705e0ffe47c07611e84a3b245" [[deps.AbstractTrees]] git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177" @@ -175,6 +175,27 @@ deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" version = "1.6.0" +[[deps.FastPower]] +git-tree-sha1 = "3d7a286bba4b6eb8904658283678be89a08da0b1" +uuid = "a4df4552-cc26-4903-aec0-212e50a0e84b" +version = "1.0.0" + + [deps.FastPower.extensions] + FastPowerEnzymeExt = "Enzyme" + FastPowerForwardDiffExt = "ForwardDiff" + FastPowerMeasurementsExt = "Measurements" + FastPowerMonteCarloMeasurementsExt = "MonteCarloMeasurements" + FastPowerReverseDiffExt = "ReverseDiff" + FastPowerTrackerExt = "Tracker" + + [deps.FastPower.weakdeps] + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" + MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + [[deps.FieldMetadata]] git-tree-sha1 = "c279c6eab9767a3f62685e5276c850512e0a1afd" uuid = "bf96fef3-21d2-5d20-8afa-0e7d4c32a885" diff --git a/Project.toml b/Project.toml index ac49bd9a5..6c896043a 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330" CFTime = "179af706-886a-5703-950a-314cd64e0468" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b" FieldMetadata = "bf96fef3-21d2-5d20-8afa-0e7d4c32a885" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" @@ -32,6 +33,7 @@ CFTime = "0.1" Dates = "<0.0.1,1" DelimitedFiles = "<0.0.1,1" Downloads = "<0.0.1,1" +FastPower = "1" FieldMetadata = "0.3" Glob = "1.3" Graphs = "1.4" diff --git a/src/Wflow.jl b/src/Wflow.jl index f4a8bd8f5..042d139f6 100644 --- a/src/Wflow.jl +++ b/src/Wflow.jl @@ -21,6 +21,7 @@ using Glob using Polyester using LoopVectorization using IfElse +using FastPower: fastpower as pow @metadata get_units "mm dt-1" String # metadata for BMI grid diff --git a/src/utils.jl b/src/utils.jl index 1f557b1cb..f0baad605 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -416,10 +416,6 @@ function det_surfacewidth(dw, riverwidth, river) return sw end -# 2.5x faster power method -"Faster method for exponentiation" -pow(x, y) = exp(y * log(x)) - "Return the sum of the array `A` at indices `inds`" function sum_at(A, inds) v = zero(eltype(A)) diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl index eb3c96763..b0fc9f05b 100644 --- a/test/horizontal_process.jl +++ b/test/horizontal_process.jl @@ -39,10 +39,10 @@ Q = zeros(n) Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, alpha, beta, DCL, dt_sec) @testset "flow rate" begin - @test sum(Q) ≈ 2.957806043289641e6 - @test Q[toposort[1]] ≈ 0.007260052312634069f0 - @test Q[toposort[n-100]] ≈ 3945.762718338739f0 - @test Q[sink] ≈ 4131.101474418251 + @test sum(Q) ≈ 2.9577875f6 + @test Q[toposort[1]] ≈ 0.0072600525f0 + @test Q[toposort[n-100]] ≈ 3945.7627f0 + @test Q[sink] ≈ 4131.0044f0 end @testset "kinematic wave subsurface flow" begin diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl index aec71a06e..c4eff55e8 100644 --- a/test/reservoir_lake.jl +++ b/test/reservoir_lake.jl @@ -53,11 +53,11 @@ end ) Wflow.update(lake, 1, 2500.0, 181, 86400.0) - @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh)[1] ≈ 19.672653848925634 - @test lake.outflow[1] ≈ 85.14292808113598 - @test lake.totaloutflow[1] ≈ 7.356348986210149e6 - @test lake.storage[1] ≈ 3.55111879238499e9 - @test lake.waterlevel[1] ≈ 19.672653848925634 + @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh)[1] ≈ 19.673122f0 + @test lake.outflow[1] ≈ 84.16451f0 + @test lake.totaloutflow[1] ≈ 7.271814f6 + @test lake.storage[1] ≈ 3.5512033f9 + @test lake.waterlevel[1] ≈ 19.673122f0 @test lake.precipitation[1] ≈ 20.0 @test lake.evaporation[1] ≈ 3.2 @test lake.actevap[1] ≈ 3.2 @@ -107,8 +107,8 @@ sh = [ Wflow.update(lake, 1, 500.0, 15, 86400.0) Wflow.update(lake, 2, 500.0, 15, 86400.0) - @test lake.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2 - @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3 + @test lake.outflow ≈ [214.77669f0, 236.83281f0] atol = 1e-2 + @test lake.totaloutflow ≈ [1.8556706f7, 2.0462356f7] atol = 1e3 @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4 @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2 lake.actevap .= 0.0 @@ -117,7 +117,7 @@ sh = [ Wflow.update(lake, 1, 500.0, 15, 86400.0) Wflow.update(lake, 2, 500.0, 15, 86400.0) @test lake.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2 - @test lake.totaloutflow ≈ [-2.2446764487487033e7, 4.3154002238515094e7] atol = 1e3 + @test lake.totaloutflow ≈ [-2.2449248f7, 4.3155704f7] atol = 1e3 @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4 @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2 @test lake.actevap ≈ [2.0, 2.0] diff --git a/test/run_sbm.jl b/test/run_sbm.jl index d50ccb752..c0b1091ce 100644 --- a/test/run_sbm.jl +++ b/test/run_sbm.jl @@ -405,7 +405,7 @@ model = Wflow.run_timestep(model) @test q[1622] ≈ 7.315676384849305f-5 @test q[43] ≈ 11.92787156357908f0 @test q[501] ≈ 3.510668846752431f0 - @test q[5808] ≈ 0.002223993845806248f0 + @test q[5808] ≈ 0.0022225645f0 h = model.lateral.river.h_av @test h[1622] ≈ 0.001987887580593841f0 @test h[43] ≈ 0.436641524481545f0 diff --git a/test/vertical_process.jl b/test/vertical_process.jl index 9ce59fb45..784ba7f3a 100644 --- a/test/vertical_process.jl +++ b/test/vertical_process.jl @@ -12,7 +12,7 @@ using Dates (4.343, 3.87, 0.387, 0.0, 3.8, 2.043), ), ) - @test Wflow.head_brooks_corey(0.25, 0.6, 0.15, 10.5, -10.0) ≈ -90.6299820833844 + @test Wflow.head_brooks_corey(0.25, 0.6, 0.15, 10.5, -10.0) ≈ -90.633095f0 @test Wflow.feddes_h3(-300.0, -600.0, 3.5, Second(86400)) ≈ -412.5 @test Wflow.feddes_h3(-300.0, -600.0, 0.5, Second(86400)) == -600.0 @test Wflow.feddes_h3(-300.0, -600.0, 6.0, Second(86400)) == -300.0 @@ -33,7 +33,7 @@ using Dates @test all( isapprox.( Wflow.unsatzone_flow_layer(43.5, 256.0, 135.0, 12.6), - (43.49983744545384, 0.00016255454615829025), + (43.499836f0, 0.00016239971f0), ), ) @test all( From fe2fc8a42e50afa168502f5493a4a7374922421a Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 1 Oct 2024 20:59:23 +0200 Subject: [PATCH 2/2] Replace squares and square roots For faster alternatives --- src/reservoir_lake.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl index 5060c842d..3ca3d1f4f 100644 --- a/src/reservoir_lake.jl +++ b/src/reservoir_lake.jl @@ -509,16 +509,16 @@ function update(lake::Lake, i, inflow, doy, timestepsecs) # outflowfunc = 3 # Calculate lake factor and SI parameter if lake.outflowfunc[i] == 3 - lakefactor = lake.area[i] / (timestepsecs * pow(lake.b[i], 0.5)) + lakefactor = lake.area[i] / (timestepsecs * sqrt(lake.b[i])) si_factor = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow # Adjust SIFactor for ResThreshold != 0 si_factor_adj = si_factor - lake.area[i] * lake.threshold[i] / timestepsecs # Calculate the new lake outflow/waterlevel/storage if si_factor_adj > 0.0 quadratic_sol_term = - -lakefactor + pow((pow(lakefactor, 2.0) + 4.0 * si_factor_adj), 0.5) + -lakefactor + sqrt(lakefactor^2 + 4.0 * si_factor_adj) if quadratic_sol_term > 0.0 - outflow = pow(0.5 * quadratic_sol_term, 2.0) + outflow = 0.5 * quadratic_sol_term^2 else outflow = 0.0 end