Hmw ct jsw/pisces

Manifest.toml

@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.2"
+julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "bd1de93833dafc8fa6de30ae101c862f7fa5877e"
+project_hash = "37e11a3bdd396c973917441db34ded05b97faa85"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -171,7 +171,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.1.1+0"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -443,20 +443,16 @@ uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e"
 version = "0.2.1+0"
 
 [[deps.KernelAbstractions]]
-deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "35ceea58aa34ad08b1ae00a52622c62d1cfb8ce2"
+deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
+git-tree-sha1 = "0fac59881e91c7233a9b0d47f4b7d9432e534f0f"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.24"
+version = "0.9.23"
 
     [deps.KernelAbstractions.extensions]
     EnzymeExt = "EnzymeCore"
-    LinearAlgebraExt = "LinearAlgebra"
-    SparseArraysExt = "SparseArrays"
 
     [deps.KernelAbstractions.weakdeps]
     EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
-    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[deps.LLVM]]
 deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
@@ -673,9 +669,9 @@ version = "1.2.0"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "12cdd648cc342e52686515811c69bc4bc452a94c"
+git-tree-sha1 = "ed415deb1d80c2a15c9b8bf0c7df04c6dec9d6c3"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.91.9"
+version = "0.91.8"
 
     [deps.Oceananigans.extensions]
     OceananigansEnzymeExt = "Enzyme"

src/Models/AdvectedPopulations/PISCES/PISCES.jl

@@ -950,5 +950,5 @@ include("zooplankton.jl")
 
 @inline sinking_tracers(::PISCES) = (:POC, :GOC, :SFe, :BFe, :PSi, :CaCO₃) # please list them here
 
-@inline chlorophyll(model) = model.tracers.Pᶜʰˡ + model.tracers.Dᶜʰˡ
+@inline chlorophyll(bgc::PISCES, model) = model.tracers.Pᶜʰˡ + model.tracers.Dᶜʰˡ
 end # module

src/Models/AdvectedPopulations/PISCES/calcite.jl

@@ -39,7 +39,7 @@ end
     wᴾ = bgc.min_quadratic_mortality_of_phytoplankton
     ηᶻ = bgc.fraction_of_calcite_not_dissolving_in_guts.Z
     ηᴹ = bgc.fraction_of_calcite_not_dissolving_in_guts.M
-    sh = shear_rate(z, zₘₓₗ)PAR₃
+    sh = shear_rate(z, zₘₓₗ)
     return rain_ratio(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, T, PAR, zₘₓₗ, bgc)*(ηᶻ*grazing_Z(P, D, POC, T, bgc)[2]*Z+ηᴹ*grazing_M(P, D, Z, POC, T, bgc)[2]*M + 0.5*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2)) #eq76
 end
 

validation/PISCES/boxPISCES.jl

@@ -18,21 +18,17 @@ using Oceananigans.Fields: FunctionField
 const year = years = 365day
 nothing #hide
 
-grid = RectilinearGrid( topology = (Flat, Flat, Flat), size = (), z = -100)
+grid = RectilinearGrid( topology = (Flat, Flat, Flat), size = (), z = 0)
 clock = Clock(time = zero(grid))
 
 # This is forced by a prescribed time-dependent photosynthetically available radiation (PAR)
 PAR_func(t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
 MLD(t) = - 100
 EU(t) = - 50
- # specify the nominal depth of the box for the PAR profile
-# Modify the PAR based on the nominal depth and exponential decay
-#PAR_func(t) = 18.0 # Modify the PAR based on the nominal depth and exponential decay
-#PAR_func(t) = 18.0
 PAR_func1(t) = PAR_func(t)/3
 PAR_func2(t) = PAR_func(t)/3
 PAR_func3(t)= PAR_func(t)/3
-
+#Each PAR component must sum to the initial PAR function at the surface, and we are only considering the 0 layer
 
 PAR₁ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
 PAR₂ = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
@@ -42,8 +38,8 @@ zₑᵤ = FunctionField{Center, Center, Center}(EU, grid; clock)
 nothing #hide
 
 # Set up the model. Here, first specify the biogeochemical model, followed by initial conditions and the start and end times
-model = BoxModel(; biogeochemistry = PISCES(; grid, light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR₁, PAR₂, PAR₃)), mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ),
-
+model = BoxModel(; biogeochemistry = PISCES(; grid, light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR₁, PAR₂, PAR₃)), 
+                   mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ),
                    clock)
 
 set!(model, P = 6.95, D = 6.95, Z = 0.695,  M = 0.695, Pᶜʰˡ = 1.671,  Dᶜʰˡ = 1.671, Pᶠᵉ =7e-6 * 1e9 / 1e6 * 6.95, Dᶠᵉ = 7e-6 * 1e9 / 1e6 * 6.95, Dˢⁱ = 1.162767, DOC = 0.0, POC = 0.0, GOC = 0.0, SFe = 7e-6 * 1e9 / 1e6 *1.256, BFe =7e-6 * 1e9 / 1e6 *1.256, NO₃ = 6.202, NH₄ = 0.25*6.202, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.29679635764590534, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = 14.0) #Using Copernicus Data (26.665, 14.)
@@ -54,12 +50,8 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, model.fields; filen
 
 prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(simulation.run_wall_time))"
 
+#NaN checker function, could be removed if confident of model stability
 function non_zero_fields!(model) 
-    #for tracer in model.fields
-      #  if tracer != model.fields[:T]
-       #     tracer = tracer[1,1,1]
-        #    getproperty(model, tracer) = max(0.0, getproperty(model, tracer))
-        #end
     @inbounds for (idx, field) in enumerate(model.fields)
         if isnan(field[1,1,1])
             throw("$(keys(model.fields)[idx]) has gone NaN")
@@ -68,9 +60,7 @@ function non_zero_fields!(model)
         end
 
     end
-    #end
     return nothing
-
 end
 
 simulation.callbacks[:progress] = Callback(prog, IterationInterval(100))
@@ -97,11 +87,13 @@ for (name, tracer) in pairs(timeseries)
     push!(axs, Axis(fig[floor(Int, idx/4), Int(idx%4)], ylabel = "$name", xlabel = "years", xticks=(0:40)))
     lines!(axs[end], times / year, tracer, linewidth = 3)
 end
+
+#Plotting the function of PAR
 push!(axs, Axis(fig[6,1], ylabel = "PAR", xlabel = "years", xticks=(0:40)))
 lines!(axs[end], (0:10day:5years) / year, x -> PAR_func(x * year), linewidth = 3)
-fi = length(timeseries.P)
 
-#println("P = $(timeseries.P[fi]), D = $(timeseries.D[fi]), Z = $(timeseries.Z[fi]),  M = $(timeseries.M[fi]), Pᶜʰˡ = $(timeseries.Pᶜʰˡ[fi]),  Dᶜʰˡ = $(timeseries.Dᶜʰˡ[fi]), Pᶠᵉ = $(timeseries.Pᶠᵉ[fi]), Dᶠᵉ = $(timeseries.Dᶠᵉ[fi]), Dˢⁱ = $(timeseries.Dˢⁱ[fi]), DOC = $(timeseries.DOC[fi]), POC = $(timeseries.POC[fi]), GOC = $(timeseries.GOC[fi]), SFe = $(timeseries.SFe[fi]), BFe = $(timeseries.BFe[fi]), PSi = $(timeseries.PSi[fi]), NO₃ = $(timeseries.NO₃[fi]), NH₄ = $(timeseries.NH₄[fi]), PO₄ = $(timeseries.PO₄[fi]), Fe = $(timeseries.Fe[fi]), Si = $(timeseries.Si[fi]), CaCO₃ = $(timeseries.CaCO₃[fi]), DIC = $(timeseries.DIC[fi]), Alk = $(timeseries.Alk[fi]), O₂ = $(timeseries.O₂[fi]), T = 14.0")
+#Below is merely a tester to check that values of things are conserved, these can be removed in the test files
+fi = length(timeseries.P)
 
 Carbon_at_start = timeseries.P[1] + timeseries.D[1] + timeseries.Z[1] + timeseries.M[1] + timeseries.DOC[1] + timeseries.POC[1] + timeseries.GOC[1] + timeseries.DIC[1] + timeseries.CaCO₃[1]
 Carbon_at_end = timeseries.P[fi] + timeseries.D[fi] + timeseries.Z[fi] + timeseries.M[fi] + timeseries.DOC[fi] + timeseries.POC[fi] + timeseries.GOC[fi] + timeseries.DIC[fi] + timeseries.CaCO₃[fi]
@@ -119,4 +111,6 @@ println("Iron at start = ", Iron_at_start, " Iron at end = ", Iron_at_end)
 println("Silicon at start = ", Silicon_at_start, " Silicon at end = ", Silicon_at_End)
 println("Phosphates at start = ", Phosphates_at_start, " Phosphates at end = ", Phosphates_at_end)
 println("Nitrogen at start = ", Nitrogen_at_start, " Nitrogen at end = ", Nitrogen_at_end)
+
+#Display the figure
 fig

validation/PISCES/columnPISCES.jl

@@ -21,9 +21,9 @@ const year = years = 365days
 nothing #hide
 
 # ## Surface PAR and turbulent vertical diffusivity based on idealised mixed layer depth 
-# Setting up idealised functions for PAR and diffusivity (details here can be ignored but these are typical of the North Atlantic)
+# Setting up idealised functions for PAR and diffusivity (details here can be ignored but these are typical of the North Atlantic), temperaeture and euphotic layer
 
-@inline PAR⁰(t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
+@inline PAR⁰(x,y,t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
 
 @inline H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)
 
@@ -33,41 +33,23 @@ nothing #hide
 
 @inline κₜ(x, y, z, t) = (1e-2 * (1 + tanh((z - MLD(x, y, z, t)) / 10)) / 2 + 1e-4)
 
-#@inline temp(x, y, z, t) = (2.4 * cos(t * 2π / year + 50days) + 10)*exp(z/50)
-@inline temp(x, y, z, t) = 14*exp(z/10)
+@inline temp(x, y, z, t) = (2.4 * cos(t * 2π / year + 50days) + 10)*exp(z/10)
 
 @inline euphotic(x, y, z, t) = - 50.0
 nothing #hide
 
-PAR_func(x, y, z, t) = PAR⁰(x,y,t)*exp(z/10) # Modify the PAR based on the nominal depth and exponential decay
-
-PAR_func1(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10)
-PAR_func2(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10)
-PAR_func3(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10)
-
-yearly_maximum_silicate = ConstantField(1) 
-dust_deposition = ConstantField(0)
-carbonate_sat_ratio = ConstantField(0)
-
-#w_GOC(z) = 30/day + (200/day - 30/day)*(max(0, abs(z)-100))/(5000)
+#The commented equation is the correct form of w_GOC, but have not figured out how to implement this
+#w_GOC(z) = 30/day + (200/day - 30/day)*(max(0, abs(z)-abs(zₘₓₗ)))/(5000)
 w_GOC = 30/day
 w_POC = 2.0/day
 grid = RectilinearGrid(size = (1, 1, 100), extent = (20meters, 20meters, 400meters))
 
 clock = Clock(; time = 0.0)
 
-PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)
-
-PAR₁ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
-PAR₂ = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
-PAR₃ = FunctionField{Center, Center, Center}(PAR_func3, grid; clock)
 zₘₓₗ = FunctionField{Center, Center, Center}(MLD, grid; clock)
 zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock)
 
 #ff = FunctionField{Nothing, Nothing, Face}(w_GOC, grid)
-# ## Grid
-# Define the grid.
-
 
 # ## Model
 # First we define the biogeochemical model including carbonate chemistry (for which we also define temperature (``T``) and salinity (``S``) fields)
@@ -77,7 +59,7 @@ zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock)
 biogeochemistry = PISCES(; grid,
                            light_attenuation_model = MultiBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR = PAR⁰), 
                            sinking_speeds = (; POC = w_POC, SFe = w_POC, GOC = w_GOC, BFe = w_GOC, PSi = w_GOC, CaCO₃ = w_GOC),
-                            mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ)
+                           mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ)
 
 CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()
 
@@ -89,17 +71,19 @@ model = NonhydrostaticModel(; grid,
                               clock,
                               closure = ScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ = κₜ),
                               biogeochemistry,
-                              boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ))
-
+                              boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ),
+                              auxiliary_fields = (; S)
+                              )
+
 
 @info "Setting initial values..."
-set!(model, P = 6.95, D = 6.95, Z = 0.695,  M = 0.695, Pᶜʰˡ = 1.671,  Dᶜʰˡ = 1.671, Pᶠᵉ =7e-6 * 1e9 / 1e6 * 6.95, Dᶠᵉ = 7e-6 * 1e9 / 1e6 * 6.95, Dˢⁱ = 1.162767, DOC = 0.0, POC = 0.0, GOC = 0.0, SFe = 7e-6 * 1e9 / 1e6 *1.256, BFe =7e-6 * 1e9 / 1e6 *1.256, NO₃ = 6.202e-2, NH₄ = 0.25*6.202e-2, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.29679635764590534, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.)
+set!(model, P = 6.95, D = 6.95, Z = 0.695,  M = 0.695, Pᶜʰˡ = 1.671,  Dᶜʰˡ = 1.671, Pᶠᵉ =7e-6 * 1e9 / 1e6 * 6.95, Dᶠᵉ = 7e-6 * 1e9 / 1e6 * 6.95, Dˢⁱ = 1.162767, DOC = 0.0, POC = 0.0, GOC = 0.0, SFe = 7e-6 * 1e9 / 1e6 *1.256, BFe =7e-6 * 1e9 / 1e6 *1.256, NO₃ = 6.202, NH₄ = 0.25*6.202, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.001, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.), Calcite is not correct, but this is to see it on the graphs
 # ## Simulation
 # Next we setup a simulation and add some callbacks that:
 # - Show the progress of the simulation
 # - Store the model and particles output
 
-simulation = Simulation(model, Δt = 90minutes, stop_time = 2years)
+simulation = Simulation(model, Δt = 50minutes, stop_time = 2years)
 
 progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n",
                                 iteration(sim),
@@ -110,6 +94,7 @@ progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time:
 simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10day)
 )
 
+#NaN Checker function. Could be removed to improve speed, if confident of model stability
 function non_zero_fields!(model) 
     @inbounds for (idx, tracer) in enumerate(model.tracers)
         for i in 1:50
@@ -179,7 +164,7 @@ carbon_export = zeros(length(times))
 using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity
 
 for (i, t) in enumerate(times)
-    air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35)
+    #air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35)
     carbon_export[i] = (POC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:POC)).w[1, 1, grid.Nz-20] +
                         GOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:GOC)).w[1, 1, grid.Nz-20]) * redfield(Val(:GOC), model.biogeochemistry)
 end
@@ -299,4 +284,8 @@ lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year,
 lines!(axfDIC, times / days, carbon_export / 1e3    * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export")
 Legend(fig[7, 2], axfDIC, framevisible = false)
 
+#Plotting a graph of Mixed Layer Depth
+axs = []
+push!(axs, Axis(fig[7,3], xlabel = "Time (days)", title = "Mixed Layer Depth (m)"))
+lines!(axs[end], (0:1day:2years)/days, x ->  MLD(0, 0, 0, x*days), linewidth = 3)
 fig