From aa10e53e636f7d988561e1d82d107d93c675feb5 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Thu, 22 Aug 2024 15:23:20 -0400 Subject: [PATCH 1/7] more iron pools --- src/Diagnostics/diagnostics_kernels.jl | 2 +- src/Diagnostics/diagnostics_struct.jl | 2 +- src/Parameters/param_default.jl | 4 ++ src/Plankton/IronEnergyMode/growth_kernels.jl | 54 +++++++++++++------ .../IronEnergyMode/plankton_generation.jl | 12 +++-- .../IronEnergyMode/plankton_growth.jl | 1 + src/Plankton/utils.jl | 8 +-- 7 files changed, 59 insertions(+), 24 deletions(-) diff --git a/src/Diagnostics/diagnostics_kernels.jl b/src/Diagnostics/diagnostics_kernels.jl index 82cc664a..fb645d58 100644 --- a/src/Diagnostics/diagnostics_kernels.jl +++ b/src/Diagnostics/diagnostics_kernels.jl @@ -20,7 +20,7 @@ function diags_spcs!(diags_sp, plank, ac, x, y, z, mode::AbstractMode, arch::Arc elseif isa(mode, MacroMolecularMode) diags = (:PS, :VDOC, :VHN4, :VNO3, :VPO4, :S_PRO, :S_DNA, :S_RNA, :resp, :ρChl, :CH, :NST, :PST, :PRO, :DNA, :RNA, :Chl) elseif isa(mode, IronEnergyMode) - diags = (:PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :BS, :VNH4, :VNO3, :VPO4, :VFe, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :Chl) + diags = (:PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :Chl) end for diag in keys(diags_sp) diff --git a/src/Diagnostics/diagnostics_struct.jl b/src/Diagnostics/diagnostics_struct.jl index 36155984..be569f30 100644 --- a/src/Diagnostics/diagnostics_struct.jl +++ b/src/Diagnostics/diagnostics_struct.jl @@ -110,7 +110,7 @@ function plank_avail_diags(mode::AbstractMode) elseif isa(mode, MacroMolecularMode) plank_avail = (:num, :graz, :mort, :dvid, :PS, :VDOC, :VNH4, :VNO3, :VPO4, :resp, :ρChl, :S_PRO, :S_DNA, :S_RNA, :exu, :CH, :NST, :PST, :PRO, :DNA, :RNA, :Chl) elseif isa(mode, IronEnergyMode) - plank_avail = (:num, :graz, :mort, :dvid, :PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :BS, :VNH4, :VNO3, :VPO4, :VFe, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :Chl) + plank_avail = (:num, :graz, :mort, :dvid, :PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :Chl) end return plank_avail end diff --git a/src/Parameters/param_default.jl b/src/Parameters/param_default.jl index 867d2ef0..9f4f80b9 100644 --- a/src/Parameters/param_default.jl +++ b/src/Parameters/param_default.jl @@ -165,6 +165,10 @@ function phyt_params_default(N::Int64, mode::IronEnergyMode) "VNH4max" => [6.9e-6], # Maximum N uptake rate (mmolN/mmolC/second) "VNO3max" => [6.9e-6], # Maximum N uptake rate (mmolN/mmolC/second) "VPO4max" => [1.2e-6], # Maximum P uptake rate (mmolP/mmolC/second) + "k_Fe_ST2PS"=> [4.6e-6], # Allocation rate of Fe from storage to PS (per second) + "k_Fe_PS2ST"=> [4.6e-6], # Allocation rate of Fe from PS to storage (per second) + "k_Fe_ST2NR"=> [4.6e-6], # Allocation rate of Fe from storage to NR (per second) + "k_Fe_NR2ST"=> [4.6e-6], # Allocation rate of Fe from NR to storage (per second) "KfePS" => [3.0e-6], # Haff-saturation coeff of iron quota for photosynthesis (mmolFe/mmolC) "KfeNR" => [2.0e-6], # Haff-saturation coeff of iron quota for NO3 reduction (mmolFe/mmolC) "KsatNH4" => [0.005], # Half-saturation coeff (mmolN/m³) diff --git a/src/Plankton/IronEnergyMode/growth_kernels.jl b/src/Plankton/IronEnergyMode/growth_kernels.jl index 1d025c9a..1e06b2e2 100644 --- a/src/Plankton/IronEnergyMode/growth_kernels.jl +++ b/src/Plankton/IronEnergyMode/growth_kernels.jl @@ -16,18 +16,11 @@ end return min(1.0f0, k/OGT_rate) end -##### intracellular iron allocation to photosynthesis -@inline function iron_alloc_PS(par, p) - f = shape_func_dec(par, p.Imax, 1.0f-2) - return f -end - ##### calculate light reaction (kJ/individual/second) -@inline function calc_PS(par, Bm, CH, Chl, qFe, p) - qFe_PS = iron_alloc_PS(par, p) * qFe - qFe_PS = qFe_PS / max(1.0f-30, Bm + CH) +@inline function calc_PS(par, Bm, CH, Chl, qFePS, p) + Qfe_ps = qFePS / max(1.0f-30, Bm + CH) αI = par * p.α * Chl / max(1.0f-30, Bm) - PS = p.PCmax * (1.0f0 - exp(-αI)) * qFe_PS / max(1.0f-30, qFe_PS + p.KfePS) * Bm + PS = p.PCmax * (1.0f0 - exp(-αI)) * Qfe_ps / max(1.0f-30, Qfe_ps + p.KfePS) * Bm return PS end @@ -161,11 +154,10 @@ function update_state_2!(plank, ΔT, arch::Architecture) end ##### nitrate reduction (mmolN/individual/second) -@inline function calc_NO3_reduction(En, qNO3, qNH4, qFe, Bm, CH, par, p, ac, ΔT) +@inline function calc_NO3_reduction(En, qNO3, qNH4, qFeNR, Bm, CH, p, ac, ΔT) reg = shape_func_dec(qNH4, p.qNH4max, 1.0f-4) - qFe_NR = max(1.0f-30, (1.0f0 - iron_alloc_PS(par, p)) * qFe) - qFe_NR = qFe_NR / max(1.0f-30, Bm + CH) - NR = p.k_nr * reg * qNO3 * qFe_NR / max(1.0f-30, qFe_NR + p.KfeNR) * ac + Qfe_NR = qFeNR / max(1.0f-30, Bm + CH) + NR = p.k_nr * reg * qNO3 * Qfe_NR / max(1.0f-30, Qfe_NR + p.KfeNR) * ac NR = min(NR, qNO3/ΔT, En/p.e_nr/ΔT) # double check qNO3 and energy are not over consumed ENR = NR * p.e_nr * ac return NR, ENR @@ -175,7 +167,7 @@ end i = @index(Global) @inbounds plank.NR[i], plank.ENR[i] = calc_NO3_reduction(plank.En[i], plank.qNO3[i], plank.qNH4[i], plank.qFe[i], plank.Bm[i], plank.CH[i], - nuts.par[i], p, plank.ac[i], ΔT) + p, plank.ac[i], ΔT) end function calc_NO3_reduc!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_NO3_reduc_kernel!(device(arch), 256, (size(plank.ac,1))) @@ -221,6 +213,34 @@ function calc_BS!(plank, p, arch::Architecture) return nothing end +##### iron allocation +@inline function iron_alloc(par, qFe, qFePS, qFeNR, p, ΔT) + reg = shape_func_dec(par, p.Imax, 1.0f-2; pow = 6.0f0) + f_ST2PS = p.k_Fe_ST2PS * reg * qFe + f_PS2ST = p.k_Fe_PS2ST * (1.0f0 - reg) * qFePS #* isless(p.Imax, par) + f_ST2NR = p.k_Fe_ST2NR * (1.0f0 - reg) * qFe + f_NR2ST = p.k_Fe_NR2ST * reg * qFeNR #* isless(par, p.Imax) + + f_ST2PS = min(f_ST2PS, qFe/ΔT) + f_PS2ST = min(f_PS2ST, qFePS/ΔT) + f_ST2NR = min(f_ST2NR, qFe/ΔT - f_ST2PS) # photosynthesis has the highest priority + f_NR2ST = min(f_NR2ST, qFeNR/ΔT) + + return f_ST2PS, f_PS2ST, f_ST2NR, f_NR2ST +end + +@kernel function calc_iron_fluxes_kernel!(plank, nuts, p, ΔT) + i = @index(Global) + @inbounds plank.ST2PS[i], plank.PS2ST[i], plank.ST2NR[i], plank.NR2ST[i] = + iron_alloc(nuts.par[i], plank.qFe[i], plank.qFePS[i], plank.qFeNR[i], + p, ΔT) +end +function calc_iron_fluxes!(plank, nuts, p, ΔT, arch::Architecture) + kernel! = calc_iron_fluxes_kernel!(device(arch), 256, (size(plank.ac,1))) + kernel!(plank, nuts, p, ΔT) + return nothing +end + ##### update C, N, P quotas, biomass, Chla @kernel function update_biomass_kernel!(plank, p, ΔT) i = @index(Global) @@ -230,6 +250,10 @@ end @inbounds plank.qP[i] -= ΔT * plank.BS[i] * p.R_PC @inbounds plank.Chl[i] += ΔT * plank.BS[i] * p.R_NC * plank.ρChl[i] @inbounds plank.age[i] += ΔT / 3600.0f0 * plank.ac[i] + @inbounds plank.qFe[i] += ΔT * (plank.PS2ST[i] - plank.ST2PS[i] + + plank.NR2ST[i] - plank.ST2NR[i]) + @inbounds plank.qFePS[i]+= ΔT * (plank.ST2PS[i] - plank.PS2ST[i]) + @inbounds plank.qFeNR[i]+= ΔT * (plank.ST2NR[i] - plank.NR2ST[i]) end function update_biomass!(plank, p, ΔT, arch::Architecture) kernel! = update_biomass_kernel!(device(arch), 256, (size(plank.ac,1))) diff --git a/src/Plankton/IronEnergyMode/plankton_generation.jl b/src/Plankton/IronEnergyMode/plankton_generation.jl index 9e14b720..3b901540 100644 --- a/src/Plankton/IronEnergyMode/plankton_generation.jl +++ b/src/Plankton/IronEnergyMode/plankton_generation.jl @@ -3,12 +3,14 @@ function construct_plankton(arch::Architecture, sp::Int, params::Dict, maxN::Int xi = zeros(Int,maxN), yi = zeros(Int,maxN), zi = zeros(Int,maxN), Sz = zeros(FT, maxN), Bm = zeros(FT, maxN), En = zeros(FT, maxN), CH = zeros(FT, maxN), qNH4 = zeros(FT, maxN), qNO3 = zeros(FT, maxN), - qFe = zeros(FT, maxN), qP = zeros(FT, maxN), Chl = zeros(FT, maxN), - gen = zeros(FT, maxN), age = zeros(FT, maxN), ac = zeros(FT, maxN), - idx = zeros(Int,maxN), + qFe = zeros(FT, maxN), qFePS= zeros(FT, maxN), qFeNR= zeros(FT, maxN), + qP = zeros(FT, maxN), Chl = zeros(FT, maxN), gen = zeros(FT, maxN), + age = zeros(FT, maxN), ac = zeros(FT, maxN), idx = zeros(Int,maxN), PS = zeros(FT, maxN), CF = zeros(FT, maxN), ECF = zeros(FT, maxN), VNH4 = zeros(FT, maxN), VNO3 = zeros(FT, maxN), VPO4 = zeros(FT, maxN), VFe = zeros(FT, maxN), ρChl = zeros(FT, maxN), + PS2ST= zeros(FT, maxN), ST2PS= zeros(FT, maxN), + NR2ST= zeros(FT, maxN), ST2NR= zeros(FT, maxN), RS = zeros(FT, maxN), ERS = zeros(FT, maxN), BS = zeros(FT, maxN), NR = zeros(FT, maxN), ENR = zeros(FT, maxN), graz = zeros(FT, maxN), mort = zeros(FT, maxN), dvid = zeros(FT, maxN) @@ -18,6 +20,7 @@ function construct_plankton(arch::Architecture, sp::Int, params::Dict, maxN::Int param_names=(:Nsuper, :Cquota, :SA, :mean, :var, :Chl2Cint, :α, :Topt, :Tmax, :Ea, :Imax, :PCmax, :k_cf, :e_cf, :VNO3max, :VNH4max, :VPO4max, + :k_Fe_ST2PS, :k_Fe_PS2ST, :k_Fe_ST2NR, :k_Fe_NR2ST, :KfePS, :KfeNR, :KsatNH4, :KsatNO3, :KsatPO4, :KSAFe, :CHmax, :qNH4max, :qNO3max, :qPmax, :qFemax, :Chl2N, :R_NC, :R_PC, :k_mtb, :k_rs, :e_rs, :k_nr, :e_nr, @@ -76,6 +79,9 @@ function generate_plankton!(plank, N::Int, g::AbstractGrid, arch::Architecture) plank.data.qP .= plank.data.qP .* (pqmax .* (plank.data.Bm .+ plank.data.CH) .- plank.data.Bm .* R_PC) # Pq plank.data.qFe .= plank.data.qFe .* (feqmax .* (plank.data.Bm .+ plank.data.CH)) # Fe + plank.data.qFePS.= plank.data.qFe .* 0.4 # Fe + plank.data.qFeNR.= plank.data.qFe .* 0.3 # Fe + plank.data.qFe .= plank.data.qFe .* 0.3 # Fe plank.data.Chl .= plank.data.Bm .* Chl2Cint # Chl mask_individuals!(plank.data, g, N, arch) diff --git a/src/Plankton/IronEnergyMode/plankton_growth.jl b/src/Plankton/IronEnergyMode/plankton_growth.jl index 837ab9de..2fce1ecc 100644 --- a/src/Plankton/IronEnergyMode/plankton_growth.jl +++ b/src/Plankton/IronEnergyMode/plankton_growth.jl @@ -10,6 +10,7 @@ function plankton_growth!(plank, nuts, rnd, p, ΔT, t, arch::Architecture) update_state_3!(plank, ΔT, arch) calc_ρChl!(plank, nuts.par, p, arch) calc_BS!(plank, p, arch) + calc_iron_fluxes!(plank, nuts, p, ΔT, arch) update_biomass!(plank, p, ΔT, arch) update_cellsize!(plank, p, arch) diff --git a/src/Plankton/utils.jl b/src/Plankton/utils.jl index 9867e435..ef935e99 100644 --- a/src/Plankton/utils.jl +++ b/src/Plankton/utils.jl @@ -82,15 +82,15 @@ function mask_individuals!(plank, g::AbstractGrid, N, arch) end ##### shape function - decrease from 1.0 to 0.0 while x increase from 0.0 to 1.0 -@inline function shape_func_dec(x, xmax, k) +@inline function shape_func_dec(x, xmax, k; pow = 4.0f0) fx = max(0.0f0, min(1.0f0, 1.0f0 - x / xmax)) - reg = fx^4.0f0 / (k + fx^4.0f0) + reg = fx^pow / (k + fx^pow) return reg end ##### shape function - increase from 0.0 to 1.0 while x increase from 0.0 to 1.0 -@inline function shape_func_inc(x, xmax, k) +@inline function shape_func_inc(x, xmax, k; pow = 4.0f0) fx = max(0.0f0, min(1.0f0, 1.0f0 - x / xmax)) - reg = fx^4.0f0 / (k + fx^4.0f0) + reg = fx^pow / (k + fx^pow) return 1.0f0 - reg end \ No newline at end of file From c77091fb9c7236cc5624a0e69b0a60537c06261f Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Thu, 22 Aug 2024 15:35:43 -0400 Subject: [PATCH 2/7] add more information when individuals esceed capacity --- src/Plankton/CarbonMode/plankton_update.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Plankton/CarbonMode/plankton_update.jl b/src/Plankton/CarbonMode/plankton_update.jl index b47897e3..604a594c 100644 --- a/src/Plankton/CarbonMode/plankton_update.jl +++ b/src/Plankton/CarbonMode/plankton_update.jl @@ -27,7 +27,7 @@ function plankton_update!(plank, nuts, rnd, p, plk, diags_spcs, ΔT, t, arch::Ar dvidnum = dot(plank.dvid, plank.ac) deactive_ind = findall(x -> x == 0.0f0, plank.ac) if dvidnum > length(deactive_ind) - throw(ArgumentError("number of individual exceeds the capacity")) + throw(ArgumentError("number of individual exceeds the capacity at timestep $t")) end divide!(plank, nuts, deactive_ind, arch) plank.idx .= 0 From b57d6c6cc1600a8380fda327ec9e6d82a52fc36b Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Thu, 22 Aug 2024 21:44:10 -0400 Subject: [PATCH 3/7] more specific error message --- src/Plankton/CarbonMode/plankton_update.jl | 2 +- src/Plankton/IronEnergyMode/plankton_update.jl | 2 +- src/Plankton/MacroMolecularMode/plankton_update.jl | 2 +- src/Plankton/QuotaMode/plankton_update.jl | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Plankton/CarbonMode/plankton_update.jl b/src/Plankton/CarbonMode/plankton_update.jl index 604a594c..6be0ac73 100644 --- a/src/Plankton/CarbonMode/plankton_update.jl +++ b/src/Plankton/CarbonMode/plankton_update.jl @@ -27,7 +27,7 @@ function plankton_update!(plank, nuts, rnd, p, plk, diags_spcs, ΔT, t, arch::Ar dvidnum = dot(plank.dvid, plank.ac) deactive_ind = findall(x -> x == 0.0f0, plank.ac) if dvidnum > length(deactive_ind) - throw(ArgumentError("number of individual exceeds the capacity at timestep $t")) + throw(ArgumentError("number of individual exceeds the capacity at timestep $(t/86400.0) days")) end divide!(plank, nuts, deactive_ind, arch) plank.idx .= 0 diff --git a/src/Plankton/IronEnergyMode/plankton_update.jl b/src/Plankton/IronEnergyMode/plankton_update.jl index 26016878..b46bbd9b 100644 --- a/src/Plankton/IronEnergyMode/plankton_update.jl +++ b/src/Plankton/IronEnergyMode/plankton_update.jl @@ -28,7 +28,7 @@ function plankton_update!(plank, nuts, rnd, p, plk, diags_spcs, ΔT, t, arch::Ar dvidnum = dot(plank.dvid, plank.ac) deactive_ind = findall(x -> x == 0.0f0, plank.ac) if dvidnum > length(deactive_ind) - throw(ArgumentError("number of individual exceeds the capacity")) + throw(ArgumentError("number of individual exceeds the capacity at timestep $(t/86400.0) days")) end divide!(plank, nuts, deactive_ind, arch) plank.idx .= 0 diff --git a/src/Plankton/MacroMolecularMode/plankton_update.jl b/src/Plankton/MacroMolecularMode/plankton_update.jl index cbf85bc9..2de2267c 100644 --- a/src/Plankton/MacroMolecularMode/plankton_update.jl +++ b/src/Plankton/MacroMolecularMode/plankton_update.jl @@ -28,7 +28,7 @@ function plankton_update!(plank, nuts, rnd, p, plk, diags_spcs, ΔT, t, arch::Ar dvidnum = dot(plank.dvid, plank.ac) deactive_ind = findall(x -> x == 0.0f0, plank.ac) if dvidnum > length(deactive_ind) - throw(ArgumentError("number of individual exceeds the capacity")) + throw(ArgumentError("number of individual exceeds the capacity at timestep $(t/86400.0) days")) end divide!(plank, nuts, deactive_ind, arch) plank.idx .= 0 diff --git a/src/Plankton/QuotaMode/plankton_update.jl b/src/Plankton/QuotaMode/plankton_update.jl index cbf85bc9..2de2267c 100644 --- a/src/Plankton/QuotaMode/plankton_update.jl +++ b/src/Plankton/QuotaMode/plankton_update.jl @@ -28,7 +28,7 @@ function plankton_update!(plank, nuts, rnd, p, plk, diags_spcs, ΔT, t, arch::Ar dvidnum = dot(plank.dvid, plank.ac) deactive_ind = findall(x -> x == 0.0f0, plank.ac) if dvidnum > length(deactive_ind) - throw(ArgumentError("number of individual exceeds the capacity")) + throw(ArgumentError("number of individual exceeds the capacity at timestep $(t/86400.0) days")) end divide!(plank, nuts, deactive_ind, arch) plank.idx .= 0 From 186c71a6fe4cb9665a101e6ba358833113417ac2 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Fri, 23 Aug 2024 14:17:04 -0400 Subject: [PATCH 4/7] bug fix --- src/Simulation/update.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Simulation/update.jl b/src/Simulation/update.jl index 96f92024..51684b92 100644 --- a/src/Simulation/update.jl +++ b/src/Simulation/update.jl @@ -31,8 +31,8 @@ function update!(sim::PlanktonSimulation; time_offset = (vels = false, PARF = fa end else for i in 1:sim.iterations - model_t_PARF = time_offset.PARF ? sim.model.t : (i-1)*sim.ΔT - model_t_temp = time_offset.temp ? sim.model.t : (i-1)*sim.ΔT + model_t_PARF = time_offset.PARF ? (i-1)*sim.ΔT : sim.model.t + model_t_temp = time_offset.temp ? (i-1)*sim.ΔT : sim.model.t t_par = floor(Int,model_t_PARF/sim.input.ΔT_PAR)+1 # starting from 1 copyto!(sim.model.timestepper.PARF, sim.input.PARF[:,:,t_par]) From 8dee2f464327e431facd0d26e001c69d2d602fa5 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Sat, 31 Aug 2024 21:25:57 -0400 Subject: [PATCH 5/7] bug fix for species-specific grazing --- src/Model/time_step.jl | 18 ++++++++++++------ src/Plankton/utils.jl | 18 ++++++++++++++---- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/src/Model/time_step.jl b/src/Model/time_step.jl index 94df1c6c..5c70ea04 100644 --- a/src/Model/time_step.jl +++ b/src/Model/time_step.jl @@ -33,8 +33,10 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) #### calculate accumulated Chla quantity (not concentration) and population find_inds!(model.individuals.phytos[sp].data, model.grid, model.arch) - acc_counts!(model.timestepper.Chl, model.timestepper.pop, - model.individuals.phytos[sp].data.Chl, model.individuals.phytos[sp].data.ac, + acc_chl!(model.timestepper.Chl, model.individuals.phytos[sp].data.Chl, + model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + acc_counts!(model.timestepper.pop, model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) end @@ -45,7 +47,7 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, model.nutrients.FeT.data, model.timestepper.par, model.timestepper.temp, model.timestepper.pop, model.arch) - + plankton_update!(model.individuals.phytos[sp].data, model.timestepper.nuts, model.timestepper.rnd, model.individuals.phytos[sp].p, model.timestepper.plk, diags.plankton[sp], ΔT, model.t, model.arch, model.mode) @@ -61,12 +63,16 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) #### calculate accumulated Chla quantity (not concentration) and population find_inds!(model.individuals.phytos[sp].data, model.grid, model.arch) - acc_counts!(model.timestepper.Chl, model.timestepper.pop, - model.individuals.phytos[sp].data.Chl, model.individuals.phytos[sp].data.ac, + acc_chl!(model.timestepper.Chl, model.individuals.phytos[sp].data.Chl, + model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + end + for sp in keys(model.individuals.phytos) + acc_counts!(model.timestepper.pop, model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) - find_NPT!(model.timestepper.nuts, model.individuals.phytos[sp].data.xi, + find_NPT!(model.timestepper.nuts, model.individuals.phytos[sp].data.xi, model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data, model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, diff --git a/src/Plankton/utils.jl b/src/Plankton/utils.jl index ef935e99..9f1f5f4e 100644 --- a/src/Plankton/utils.jl +++ b/src/Plankton/utils.jl @@ -39,15 +39,25 @@ function find_NPT!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, temp, pop, a return nothing end -##### calculate Chla and individual counts based on the status of plankton individuals -@kernel function acc_counts_kernel!(ctsChl, ctspop, Chl, ac, x, y, z) +##### calculate Chla based on the status of plankton individuals +@kernel function acc_chl_kernel!(ctsChl, Chl, ac, x, y, z) i = @index(Global) @inbounds KernelAbstractions.@atomic ctsChl[x[i], y[i], z[i]] += Chl[i] * ac[i] +end +function acc_chl!(ctsChl, Chl, ac, x, y, z, arch) + kernel! = acc_chl_kernel!(device(arch), 256, (size(ac,1))) + kernel!(ctsChl, Chl, ac, x, y, z) + return nothing +end + +##### calculate individual counts based on the status of plankton individuals +@kernel function acc_counts_kernel!(ctspop, ac, x, y, z) + i = @index(Global) @inbounds KernelAbstractions.@atomic ctspop[x[i], y[i], z[i]] += ac[i] end -function acc_counts!(ctsChl, ctspop, Chl, ac, x, y, z, arch) +function acc_counts!(ctspop, ac, x, y, z, arch) kernel! = acc_counts_kernel!(device(arch), 256, (size(ac,1))) - kernel!(ctsChl, ctspop, Chl, ac, x, y, z) + kernel!(ctspop, ac, x, y, z) return nothing end From cb518e44546594c612593a6d6f628aa4777a43f0 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Sat, 31 Aug 2024 22:20:04 -0400 Subject: [PATCH 6/7] more detailed iron allocation --- src/Model/time_step.jl | 48 +++++++++++-------- src/Model/timestepper.jl | 8 ++-- src/Parameters/param_default.jl | 10 ++-- src/Plankton/IronEnergyMode/growth_kernels.jl | 16 +++---- src/Plankton/Plankton.jl | 2 +- src/Plankton/utils.jl | 7 +-- 6 files changed, 50 insertions(+), 41 deletions(-) diff --git a/src/Model/time_step.jl b/src/Model/time_step.jl index 5c70ea04..5d2f51e4 100644 --- a/src/Model/time_step.jl +++ b/src/Model/time_step.jl @@ -34,23 +34,28 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) #### calculate accumulated Chla quantity (not concentration) and population find_inds!(model.individuals.phytos[sp].data, model.grid, model.arch) acc_chl!(model.timestepper.Chl, model.individuals.phytos[sp].data.Chl, - model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, - model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) acc_counts!(model.timestepper.pop, model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + ##### calculate PAR + for ki in 1:model.grid.Nz + calc_par!(model.timestepper.par, model.arch, model.timestepper.Chl, model.timestepper.PARF, + model.grid, model.bgc_params["kc"], model.bgc_params["kw"], ki) + end end for sp in keys(model.individuals.phytos) find_NPT!(model.timestepper.nuts, model.individuals.phytos[sp].data.xi, - model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, - model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data, - model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, - model.nutrients.FeT.data, model.timestepper.par, model.timestepper.temp, - model.timestepper.pop, model.arch) + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, + model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data, + model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, + model.nutrients.FeT.data, model.timestepper.par, model.timestepper.par₀, + model.timestepper.temp, model.timestepper.pop, model.arch) plankton_update!(model.individuals.phytos[sp].data, model.timestepper.nuts, - model.timestepper.rnd, model.individuals.phytos[sp].p, - model.timestepper.plk, diags.plankton[sp], ΔT, model.t, model.arch, model.mode) + model.timestepper.rnd, model.individuals.phytos[sp].p, + model.timestepper.plk, diags.plankton[sp], ΔT, model.t, model.arch, model.mode) end else # model.bgc_params["shared_graz"] ≠ 1.0 - species-specific grazing for sp in keys(model.individuals.phytos) @@ -64,8 +69,13 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) #### calculate accumulated Chla quantity (not concentration) and population find_inds!(model.individuals.phytos[sp].data, model.grid, model.arch) acc_chl!(model.timestepper.Chl, model.individuals.phytos[sp].data.Chl, - model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, - model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + model.individuals.phytos[sp].data.ac, model.individuals.phytos[sp].data.xi, + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, model.arch) + ##### calculate PAR + for ki in 1:model.grid.Nz + calc_par!(model.timestepper.par, model.arch, model.timestepper.Chl, model.timestepper.PARF, + model.grid, model.bgc_params["kc"], model.bgc_params["kw"], ki) + end end for sp in keys(model.individuals.phytos) acc_counts!(model.timestepper.pop, model.individuals.phytos[sp].data.ac, @@ -73,11 +83,11 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) model.individuals.phytos[sp].data.zi, model.arch) find_NPT!(model.timestepper.nuts, model.individuals.phytos[sp].data.xi, - model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, - model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data, - model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, - model.nutrients.FeT.data, model.timestepper.par, model.timestepper.temp, - model.timestepper.pop, model.arch) + model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi, + model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data, + model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data, + model.nutrients.FeT.data, model.timestepper.par, model.timestepper.par₀, + model.timestepper.temp, model.timestepper.pop, model.arch) plankton_update!(model.individuals.phytos[sp].data, model.timestepper.nuts, model.timestepper.rnd, model.individuals.phytos[sp].p, @@ -86,11 +96,6 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) end end - ##### calculate PAR - for ki in 1:model.grid.Nz - calc_par!(model.timestepper.par, model.arch, model.timestepper.Chl, model.timestepper.PARF, - model.grid, model.bgc_params["kc"], model.bgc_params["kw"], ki) - end ##### diagnostics for nutrients @inbounds diags.tracer.PAR .+= model.timestepper.par @inbounds diags.tracer.T .+= model.timestepper.temp @@ -106,6 +111,7 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics) @inbounds model.timestepper.vel₀.u.data .= model.timestepper.vel₁.u.data @inbounds model.timestepper.vel₀.v.data .= model.timestepper.vel₁.v.data @inbounds model.timestepper.vel₀.w.data .= model.timestepper.vel₁.w.data + @inbounds model.timestepper.par₀ .= model.timestepper.par return nothing end diff --git a/src/Model/timestepper.jl b/src/Model/timestepper.jl index e2784792..006775d7 100644 --- a/src/Model/timestepper.jl +++ b/src/Model/timestepper.jl @@ -8,6 +8,7 @@ mutable struct timestepper temp::AbstractArray # a (Cu)Array to store temperature field of each timestep plk::NamedTuple # a NamedTuple same as nutrients to store interactions with individuals par::AbstractArray # a (Cu)Array to store PAR field of each timestep + par₀::AbstractArray # a (Cu)Array to store PAR field of the previous timestep Chl::AbstractArray # a (Cu)Array to store Chl field of each timestep pop::AbstractArray # a (Cu)Array to store population field of each timestep rnd::AbstractArray # a StructArray of random numbers for plankton diffusion or grazing, mortality and division. @@ -25,6 +26,7 @@ function timestepper(arch::Architecture, FT::DataType, g::AbstractGrid, maxN) plk = nutrients_init(arch, g, FT) par = zeros(FT, g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch) + par₀= zeros(FT, g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch) Chl = zeros(FT, g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch) pop = zeros(FT, g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch) @@ -42,11 +44,11 @@ function timestepper(arch::Architecture, FT::DataType, g::AbstractGrid, maxN) nuts = StructArray(NH4 = zeros(FT, maxN), NO3 = zeros(FT, maxN), PO4 = zeros(FT, maxN), DOC = zeros(FT, maxN), FeT = zeros(FT, maxN), par = zeros(FT, maxN), - T = zeros(FT, maxN), pop = zeros(FT, maxN), idc = zeros(FT, maxN), - idc_int = zeros(Int, maxN)) + T = zeros(FT, maxN), pop = zeros(FT, maxN), dpar= zeros(FT, maxN), + idc = zeros(FT, maxN), idc_int = zeros(Int, maxN)) nuts_d = replace_storage(array_type(arch), nuts) - ts = timestepper(Gcs, nut_temp, vel₀, vel½, vel₁, PARF, temp, plk, par, Chl, pop, rnd_d, velos_d, nuts_d) + ts = timestepper(Gcs, nut_temp, vel₀, vel½, vel₁, PARF, temp, plk, par, par₀, Chl, pop, rnd_d, velos_d, nuts_d) return ts end \ No newline at end of file diff --git a/src/Parameters/param_default.jl b/src/Parameters/param_default.jl index 9f4f80b9..6224f0a4 100644 --- a/src/Parameters/param_default.jl +++ b/src/Parameters/param_default.jl @@ -158,17 +158,17 @@ function phyt_params_default(N::Int64, mode::IronEnergyMode) "Topt" => [27.0], # Optimal temperature for growth (C) "Tmax" => [30.0], # Maximal temperature for growth (C) "Ea" => [5.3e4], # Free energy - "Imax" => [2500.0], # Light regulation of iron allocation (μmol photon/m²/second) + "Imax" => [3000.0], # Light regulation of iron allocation (μmol photon/m²/second) "PCmax" => [2.43e-5], # Maximum primary production rate (kJ/mmolC/second) "k_cf" => [1.0e-5], # Carbon fixation rate (per second) "e_cf" => [0.59], # Energy consumption rate of carbon fixation (kJ/mmolC) "VNH4max" => [6.9e-6], # Maximum N uptake rate (mmolN/mmolC/second) "VNO3max" => [6.9e-6], # Maximum N uptake rate (mmolN/mmolC/second) "VPO4max" => [1.2e-6], # Maximum P uptake rate (mmolP/mmolC/second) - "k_Fe_ST2PS"=> [4.6e-6], # Allocation rate of Fe from storage to PS (per second) - "k_Fe_PS2ST"=> [4.6e-6], # Allocation rate of Fe from PS to storage (per second) - "k_Fe_ST2NR"=> [4.6e-6], # Allocation rate of Fe from storage to NR (per second) - "k_Fe_NR2ST"=> [4.6e-6], # Allocation rate of Fe from NR to storage (per second) + "k_Fe_ST2PS"=> [2.4e-5], # Allocation rate of Fe from storage to PS (per second) + "k_Fe_PS2ST"=> [1.2e-6], # Allocation rate of Fe from PS to storage (per second) + "k_Fe_ST2NR"=> [1.2e-5], # Allocation rate of Fe from storage to NR (per second) + "k_Fe_NR2ST"=> [1.2e-5], # Allocation rate of Fe from NR to storage (per second) "KfePS" => [3.0e-6], # Haff-saturation coeff of iron quota for photosynthesis (mmolFe/mmolC) "KfeNR" => [2.0e-6], # Haff-saturation coeff of iron quota for NO3 reduction (mmolFe/mmolC) "KsatNH4" => [0.005], # Half-saturation coeff (mmolN/m³) diff --git a/src/Plankton/IronEnergyMode/growth_kernels.jl b/src/Plankton/IronEnergyMode/growth_kernels.jl index 1e06b2e2..89a80604 100644 --- a/src/Plankton/IronEnergyMode/growth_kernels.jl +++ b/src/Plankton/IronEnergyMode/growth_kernels.jl @@ -214,12 +214,12 @@ function calc_BS!(plank, p, arch::Architecture) end ##### iron allocation -@inline function iron_alloc(par, qFe, qFePS, qFeNR, p, ΔT) - reg = shape_func_dec(par, p.Imax, 1.0f-2; pow = 6.0f0) - f_ST2PS = p.k_Fe_ST2PS * reg * qFe - f_PS2ST = p.k_Fe_PS2ST * (1.0f0 - reg) * qFePS #* isless(p.Imax, par) - f_ST2NR = p.k_Fe_ST2NR * (1.0f0 - reg) * qFe - f_NR2ST = p.k_Fe_NR2ST * reg * qFeNR #* isless(par, p.Imax) +@inline function iron_alloc(par, dpar, qFe, qFePS, qFeNR, p, ΔT) + reg = shape_func_dec(par, p.Imax, 1.0f-2; pow = 4.0f0) + f_ST2PS = p.k_Fe_ST2PS * reg * qFe * isless(0.0f0, dpar) + f_PS2ST = p.k_Fe_PS2ST * qFePS * (reg * isless(dpar, 0.0f0) + isequal(0.0f0, par)) + f_ST2NR = p.k_Fe_ST2NR * qFe * isequal(0.0f0, par) + f_NR2ST = p.k_Fe_NR2ST * qFeNR * isless(0.0f0, par) f_ST2PS = min(f_ST2PS, qFe/ΔT) f_PS2ST = min(f_PS2ST, qFePS/ΔT) @@ -232,8 +232,8 @@ end @kernel function calc_iron_fluxes_kernel!(plank, nuts, p, ΔT) i = @index(Global) @inbounds plank.ST2PS[i], plank.PS2ST[i], plank.ST2NR[i], plank.NR2ST[i] = - iron_alloc(nuts.par[i], plank.qFe[i], plank.qFePS[i], plank.qFeNR[i], - p, ΔT) + iron_alloc(nuts.par[i], nuts.dpar[i], plank.qFe[i], + plank.qFePS[i], plank.qFeNR[i], p, ΔT) end function calc_iron_fluxes!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_iron_fluxes_kernel!(device(arch), 256, (size(plank.ac,1))) diff --git a/src/Plankton/Plankton.jl b/src/Plankton/Plankton.jl index 7119dd16..21c3dc4e 100644 --- a/src/Plankton/Plankton.jl +++ b/src/Plankton/Plankton.jl @@ -4,7 +4,7 @@ export plankton_advection! export plankton_diffusion! export plankton_update! export generate_individuals, individuals -export find_inds!, find_NPT!, acc_counts!, calc_par! +export find_inds!, find_NPT!, acc_counts!, acc_chl!, calc_par! using StructArrays using Random diff --git a/src/Plankton/utils.jl b/src/Plankton/utils.jl index 9f1f5f4e..bc880004 100644 --- a/src/Plankton/utils.jl +++ b/src/Plankton/utils.jl @@ -22,7 +22,7 @@ function find_inds!(plank, g::AbstractGrid, arch::Architecture) return nothing end -@kernel function find_NPT_kernel!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, temp, pop) +@kernel function find_NPT_kernel!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, par₀, temp, pop) i = @index(Global) @inbounds nuts.NH4[i] = max(0.0f0, NH4[x[i], y[i], z[i]]) * ac[i] @inbounds nuts.NO3[i] = max(0.0f0, NO3[x[i], y[i], z[i]]) * ac[i] @@ -30,12 +30,13 @@ end @inbounds nuts.DOC[i] = max(0.0f0, DOC[x[i], y[i], z[i]]) * ac[i] @inbounds nuts.FeT[i] = max(0.0f0, FeT[x[i], y[i], z[i]]) * ac[i] @inbounds nuts.par[i] = par[x[i], y[i], z[i]] * ac[i] + @inbounds nuts.dpar[i]= (par[x[i], y[i], z[i]] - par₀[x[i], y[i], z[i]]) * ac[i] @inbounds nuts.T[i] =temp[x[i], y[i], z[i]] * ac[i] @inbounds nuts.pop[i] = pop[x[i], y[i], z[i]] * ac[i] end -function find_NPT!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, temp, pop, arch::Architecture) +function find_NPT!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, par₀, temp, pop, arch::Architecture) kernel! = find_NPT_kernel!(device(arch), 256, (size(ac,1))) - kernel!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, temp, pop) + kernel!(nuts, x, y, z, ac, NH4, NO3, PO4, DOC, FeT, par, par₀, temp, pop) return nothing end From 6a955e36b93a7becc938737057252cc8815991a4 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Sat, 31 Aug 2024 22:24:53 -0400 Subject: [PATCH 7/7] fix test --- test/example_3D_test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/example_3D_test.jl b/test/example_3D_test.jl index 770a294d..a05ca216 100644 --- a/test/example_3D_test.jl +++ b/test/example_3D_test.jl @@ -47,5 +47,5 @@ TCt = tot_mass(model.nutrients.DIC.data, grid) + TCt = TCt + sum(model.individuals.phytos.sp1.data.Bm) @testset "PlanktonIndividuals 3D tests:" begin - @test isapprox(TC,TCt; atol=1e2) + @test isapprox(TC,TCt; atol=1e3) end