Use FastPower.jl for faster power function

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"

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"

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

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

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))

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

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]

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

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(