diff --git a/examples/BEPS_compare_with_C.qmd b/examples/BEPS_compare_with_C.qmd new file mode 100644 index 0000000..73d1622 --- /dev/null +++ b/examples/BEPS_compare_with_C.qmd @@ -0,0 +1,61 @@ +# BEPS + +```{julia} +# import RTableTools: fread +# d = fread("examples/input/p1_meteo") +``` + +```{julia} +using BEPS +d = deserialize("data/p1_meteo") +lai = readdlm("examples/input/p1_lai.txt")[:] + +par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85, + soil_type=8, Tsoil=2.2, + soilwater=0.4115, snowdepth=0.0) +``` + +```{julia} +@time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); +@time df_c, df_ET_c = besp_main(d, lai, par; version="c"); + +sum(df_jl) +sum(df_c) + +df_diff = abs.(df_jl .- df_c) +df_diff_perc = abs.(df_jl .- df_c) ./ df_c .* 100 +# df_diff = abs.(df_ET_jl .- df_ET_c) +# fwrite(df_diff, "a-3.csv") +max(df_diff) +max(df_diff_perc) # SH, 1.48%的误差 +``` + +# 绘图展示结果 +```{julia} +using Plots +using Ipaper + +function plot_var(var=:SH) + n = size(df_jl, 1) + inds = 1:n + + y_jl = df_jl[inds, var] + y_c = df_c[inds, var] + diff_perc = (y_jl .- y_c) ./ y_c .* 100 + + plot(diff_perc; title="$var", label="bias (%)") # +end + +vars = names(df_jl)[[1:4; 8; 12:16]] +ps = map(plot_var, vars) +plot(ps..., layout=(3, 4), size=(1400, 800)) +# write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false) +``` + +```{julia} +# using ProfileView +# @time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); +# @profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); +@profview_allocs df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); +@profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); +``` diff --git a/examples/example_01.qmd b/examples/example_01.qmd index 91bb4cd..4564f2d 100644 --- a/examples/example_01.qmd +++ b/examples/example_01.qmd @@ -24,127 +24,32 @@ sum(df_c) df_diff = abs.(df_jl .- df_c) df_diff_perc = abs.(df_jl .- df_c) ./ df_c .* 100 -# df_diff = abs.(df_ET_jl .- df_ET_c) -# fwrite(df_diff, "a-3.csv") + max(df_diff) max(df_diff_perc) # SH, 1.48%的误差 ``` # 绘图展示结果 + ```{julia} using Plots using Ipaper +gr(framestyle=:box) function plot_var(var=:SH) n = size(df_jl, 1) + n = 24*7 inds = 1:n y_jl = df_jl[inds, var] - y_c = df_c[inds, var] - diff_perc = (y_jl .- y_c) ./ y_c .* 100 - - plot(diff_perc; title="$var", label="bias (%)") # + # y_c = df_c[inds, var] + # diff_perc = (y_jl .- y_c) ./ y_c .* 100 + plot(y_jl; title="$var", label=nothing) # end vars = names(df_jl)[[1:4; 8; 12:16]] ps = map(plot_var, vars) -# p1 = plot(df_jl.SH[inds], label="julia", title="SH") -# plot!(df_c.SH[inds] + 0, label="c") -# p2 = plot(df_diff_perc.SH[inds], label="Julia error (%)") -# title = "bias (%): (julia - c)/c" -plot(ps..., layout=(3, 4), size=(1400, 800)) -write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false) -``` - -```{julia} -# using ProfileView -# @time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); -# @profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); -@profview_allocs df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); -@profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia"); -``` - - -```{julia} -## C版本的结果 -# :gpp_o_sunlit => 41704.09054753295 -# :gpp_u_sunlit => 285.1932406544977 -# :gpp_o_shaded => 12720.54113242912 -# :gpp_u_shaded => 135.13475998737647 -# :plant_resp => 0.0 -# :npp_o => 0.0 -# :npp_u => 0.0 -# :GPP => 2369.3022582020903 -# :NPP => 0.0 -# :NEP => 0.0 -# :soil_resp => 0.0 -# :Net_Rad => 1.299620168997007e6 -# :SH => 232161.16420363513 -# :LH => 811896.0668881282 -# :Trans => 444.02437414775443 -# :Evap => 748.8880942489606 -## sum -:gpp_o_sunlit => 0.0002343150940907193 -:gpp_u_sunlit => 0.05115824940432174 -:gpp_o_shaded => 1.4900267617763002e-5 -:gpp_u_shaded => 1.0758724368703184e-7 -:plant_resp => 0.0 -:npp_o => 0.0 -:npp_u => 0.0 -:GPP => 0.0022203403965193147 -:NPP => 0.0 -:NEP => 0.0 -:soil_resp => 0.0 -:Net_Rad => 0.003708570688994328 -:SH => 0.12619010825226695 -:LH => 2.5288028670446403 -:Trans => 0.0035486520015893166 -:Evap => 0.0009524622286666777 - -## max - :gpp_o_sunlit => 0.00021630367300851105 - :gpp_u_sunlit => 0.0005586954093889651 - :gpp_o_shaded => 9.06128597910616e-6 - :gpp_u_shaded => 1.7276561885964936e-8 - :plant_resp => 0.0 - :npp_o => 0.0 - :npp_u => 0.0 - :GPP => 2.4137032993143404e-5 - :NPP => 0.0 - :NEP => 0.0 - :soil_resp => 0.0 - :Net_Rad => 0.002326329983901587 - :SH => 0.029240292020906722 - :LH => 0.07646965953554741 - :Trans => 1.6816338803168907e-5 - :Evap => 0.00012368929687477104 -``` - - -# Test about - -```{julia} -canopyh_o = 2.0 -canopyh_u = 0.2 -height_wind_sp = 2.0 -clumping = 0.8 -Ta = 20.0 -wind_sp = 2. -GH_o = 100.0 -pai_o = 4.0 -pai_u = 2.0 - -ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u = - aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o, - pai_o, pai_u) -r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u) - -ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u = - aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o, - pai_o, pai_u) -r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u) - -@show r1 -@show r2; +plot(ps..., layout=(3, 4), size=(1400, 800)) +# write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false) ``` diff --git a/src/BEPS/aerodynamic_conductance.jl b/src/BEPS/aerodynamic_conductance.jl index f653d3a..829663b 100644 --- a/src/BEPS/aerodynamic_conductance.jl +++ b/src/BEPS/aerodynamic_conductance.jl @@ -5,8 +5,20 @@ # Examples ```julia -ra_o, G_o_a, G_o_b, G_u_a, G_u_b, ra_g = aerodynamic_conductance_jl( - canopy_height_o, canopy_height_u, z_wind, clumping, temp_air, wind_sp, SH_o_p, lai_o, lai_u) +canopyh_o = 2.0 +canopyh_u = 0.2 +height_wind_sp = 2.0 +clumping = 0.8 +Ta = 20.0 +wind_sp = 2. +GH_o = 100.0 +lai_o = 4.0 +lai_u = 2.0 + +ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u = + aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o, + lai_o, lai_u) +r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u) ``` """ function aerodynamic_conductance_jl(canopy_height_o::FT, canopy_height_u::FT, @@ -91,7 +103,7 @@ end function windProfile_factor(canopy_height_u, canopy_height_o, gamma, k=1.0) - exp(-gamma * (1 - canopy_height_u * k / canopy_height_o)) + exp(gamma * (1 - canopy_height_u * k / canopy_height_o)) end function cal_Nu(u::FT, nu_lower::FT)::FT diff --git a/src/beps_inter_prg.jl b/src/beps_inter_prg.jl index dd3419e..0252fa6 100644 --- a/src/beps_inter_prg.jl +++ b/src/beps_inter_prg.jl @@ -36,10 +36,12 @@ function inter_prg_jl( Tc_old = Leaf() Tc_new = Leaf() Gs_old = Leaf() - Gc = Leaf() # total conductance for CO2 from the intercellular space of the leaves to the reference height above the canopy - Gh = Leaf() # total conductance for heat transfer from the leaf surface to the reference height above the canopy - Gw = Leaf() # total conductance for water from the intercellular space of the leaves to the reference height above the canopy - Gww = Leaf() # total conductance for water from the surface of the leaves to the reference height above the canopy + + # to the reference height above the canopy + Gc = Leaf() # total conductance for CO2 from the intercellular space of the leaves + Gh = Leaf() # total conductance for heat transfer from the leaf surface + Gw = Leaf() # total conductance for water from the intercellular space of the leaves + Gww = Leaf() # total conductance for water from the surface of the leaves Gs_new = Leaf() Ac = Leaf() @@ -71,7 +73,6 @@ function inter_prg_jl( Kn = 0.3 #0.713/2.4 G_theta = 0.5 gamma = 0.066 - # double K,Vcmax0,Vcmax_sunlit,Vcmax_shaded,expr1,expr2,expr3; # double slope_Vcmax_N, leaf_N,Jmax_sunlit,Jmax_shaded; @@ -435,6 +436,5 @@ function inter_prg_jl( # total GPP . gC/m2/step mid_res.GPP = (GPP.o_sunlit + GPP.o_shaded + GPP.u_sunlit + GPP.u_shaded) * 12 * step * 0.000001 - # return nothing end diff --git a/examples/debug_Rn.qmd b/test/debug-Rn.qmd similarity index 100% rename from examples/debug_Rn.qmd rename to test/debug-Rn.qmd diff --git a/test/test-photo.qmd b/test/debug-photo.qmd similarity index 100% rename from test/test-photo.qmd rename to test/debug-photo.qmd diff --git a/test/runtests.jl b/test/runtests.jl index b51791f..5cf2e30 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using BEPS: path_proj -include("test-readparam.jl") +include("test-clang.jl") include("test-soil.jl") diff --git a/test/test-beps_main.jl b/test/test-beps_main.jl deleted file mode 100644 index 9357cc3..0000000 --- a/test/test-beps_main.jl +++ /dev/null @@ -1,47 +0,0 @@ -using BEPS -using BEPS.beps_c - -# using UnPack -# using DelimitedFiles: readdlm -# include("main_pkgs.jl") - -begin - f = "examples/input/p1_meteo.txt" - d = fread(f) - lai = readdlm("examples/input/p1_lai.txt")[:] - - par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85, - soil_type=8, Tsoil=2.2, soilwater=0.4115, snowdepth=0) - # @unpack lon, lat, landcover, clumping, soil_type, Tsoil, soilwater, snowdepth = par - - # 2.1 990 3.4124 0.6607 0.0177 0.97491 0.0213 0.0531 0.9715 12.5484 8.8816 - # LAI_yr, ann_NPP, ccd, cssd, csmd, cfsd, cfmd, csm, cm, cs, cp - - # 120.5 30.5 25 0.85 8 2.2 0.4115 0.0 - # long., lat., LC, CI, soiltxt, Tsoil, soilwater,snow-dp - "" -end - -@time df_out = besp_main(d, lai, par); -sum(df_out) - -r1 = r2 = Ref(1.0) - -# :gpp_o_sunlit => 41704.09054753295 -# :gpp_u_sunlit => 285.1932406544977 -# :gpp_o_shaded => 12720.54113242912 -# :gpp_u_shaded => 135.13475998737647 -# :plant_resp => 0.0 -# :npp_o => 0.0 -# :npp_u => 0.0 -# :GPP => 2369.3022582020903 -# :NPP => 0.0 -# :NEP => 0.0 -# :soil_resp => 0.0 -# :Net_Rad => 1.299620168997007e6 -# :SH => 232161.16420363513 -# :LH => 811896.0668881282 -# :Trans => 444.02437414775443 -# :Evap => 748.8880942489606 - -# fwrite(df_out, "p1_out.csv") diff --git a/test/test-clang.jl b/test/test-clang.jl new file mode 100644 index 0000000..1a7fb09 --- /dev/null +++ b/test/test-clang.jl @@ -0,0 +1,46 @@ +using Test +using BEPS + +@testset "readparam" begin + lcs = [1, 2, 6, 9, 13, 40, -1] + for lc = lcs + par1 = readparam(lc) + par2 = clang.readparam(lc) + @test maximum(abs.(par1 .- par2)) == 0 + end +end + +@testset "readcoef" begin + lcs = [1, 6, 13, -1] + stxts = 1:11 + + for lc = lcs, stxt = stxts + coef1 = readcoef(lc, stxt) + coef2 = clang.readcoef(lc, stxt) + @test maximum(abs.(coef1 .- coef2)) == 0 + end +end + +@testset "aerodynamic_conductance" begin + canopyh_o = 2.0 + canopyh_u = 0.2 + height_wind_sp = 2.0 + clumping = 0.8 + Ta = 20.0 + wind_sp = 2.0 + GH_o = 100.0 + pai_o = 4.0 + pai_u = 2.0 + + ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u = + aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o, + pai_o, pai_u) + r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u) + + ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u = + clang.aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o, + pai_o, pai_u) + r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u) + + @test maximum(abs.(values(r1) .- values(r2))) <= 1e-8 +end diff --git a/test/test-readparam.jl b/test/test-readparam.jl index 6242450..098be73 100644 --- a/test/test-readparam.jl +++ b/test/test-readparam.jl @@ -2,24 +2,3 @@ using Test using BEPS # import BEPS: readparam, readcoef -@testset "readparam" begin - lcs = [1, 2, 6, 9, 13, 40, -1] - - for lc = lcs - par1 = readparam(lc) - par2 = clang.readparam(lc) - @test maximum(abs.(par1 .- par2)) == 0 - end -end - - -@testset "readparam" begin - lcs = [1, 6, 13, -1] - stxts = 1:11 - - for lc = lcs, stxt = stxts - coef1 = readcoef(lc, stxt) - coef2 = clang.readcoef(lc, stxt) - @test maximum(abs.(coef1 .- coef2)) == 0 - end -end