Merge pull request #170 from OceanBioME/jsw/optional-naning-for-negat…

…ivity-protection

Added `invalid_fill_value` to `ScaleNegativeTracers`

Project.toml

@@ -1,7 +1,7 @@
 name = "OceanBioME"
 uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
 authors = ["Jago Strong-Wright <[email protected]> and contributors"]
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/index.md

@@ -25,11 +25,15 @@ julia> Pkg.add("OceanBioME")
 
 ## Running your first model
 
-As a simple example lets run a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model in a two-dimensional simulation of a buoyancy front. This example requires Oceananigans, so we install that first:
+As a simple example lets run a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model in a two-dimensional simulation of a buoyancy front. This example requires Oceananigans, so we install that first via:
 
-```@example quickstart
+```julia
 using Pkg; Pkg.add("Oceananigans")
+```
 
+and then:
+
+```@example quickstart
 using OceanBioME, Oceananigans
 using Oceananigans.Units
 

src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl

@@ -303,6 +303,7 @@ function LOBSTER(; grid,
                    open_bottom::Bool = true,
 
                    scale_negatives = false,
+                   invalid_fill_value = NaN,
 
                    particles::P = nothing,
                    modifiers::M = nothing) where {FT, LA, S, P, M}
@@ -350,7 +351,7 @@ function LOBSTER(; grid,
                                          sinking_velocities)
 
     if scale_negatives
-        scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid)
+        scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid; invalid_fill_value)
         if isnothing(modifiers)
             modifiers = scaler
         elseif modifiers isa Tuple

src/Models/AdvectedPopulations/NPZD.jl

@@ -179,6 +179,7 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid,
                                                     open_bottom::Bool = true,
 
                                                     scale_negatives = false,
+                                                    invalid_fill_value = NaN,
 
                                                     particles::P = nothing,
                                                     modifiers::M = nothing) where {FT, LA, S, P, M}
@@ -200,7 +201,7 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid,
                                                  sinking_velocities)
 
     if scale_negatives
-        scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid)
+        scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid; invalid_fill_value)
         modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler)
     end
 

src/OceanBioME.jl

@@ -126,6 +126,7 @@ update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modi
 
 function update_biogeochemical_state!(bgc::Biogeochemistry, model)
     update_biogeochemical_state!(model, bgc.modifiers)
+    synchronize(device(architecture(model)))
     update_biogeochemical_state!(model, bgc.light_attenuation)
 end
 

src/Utils/negative_tracers.jl

@@ -33,21 +33,23 @@ end
 ##### Infastructure to rescale negative values
 #####
 
-struct ScaleNegativeTracers{FA, SA, W}
-    tracers :: FA
-    scalefactors :: SA
-    warn :: W
-
-    ScaleNegativeTracers(tracers::FA, scalefactors::SA, warn::W) where {FA, SA, W} =
-        warn ? error("Warning not currently implemented") : new{FA, SA, W}(tracers, scalefactors, warn)
+struct ScaleNegativeTracers{FA, SA, FV, W}
+           tracers :: FA
+      scalefactors :: SA
+invalid_fill_value :: FV
+              warn :: W
+
+    ScaleNegativeTracers(tracers::FA, scalefactors::SA, invalid_fill_value::FV, warn::W) where {FA, SA, W, FV} =
+        warn ? error("Warning not currently implemented") : new{FA, SA, FV, W}(tracers, scalefactors, invalid_fill_value, warn)
 end
 
 adapt_structure(to, snt::ScaleNegativeTracers) = ScaleNegativeTracers(adapt(to, snt.tracers),
                                                                       adapt(to, snt.scalefactors),
+                                                                      adapt(to, snt.invalid_fill_value),
                                                                       adapt(to, snt.warn))
 
 """
-    ScaleNegativeTracers(; tracers, scalefactors = ones(length(tracers)), warn = false)
+    ScaleNegativeTracers(; tracers, scalefactors = ones(length(tracers)), warn = false, invalid_fill_value = NaN)
 
 Constructs a modifier to scale `tracers` so that none are negative. Use like:
 ```julia
@@ -60,14 +62,27 @@ github](https://github.com/OceanBioME/OceanBioME.jl/discussions/48).
 
 Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative.
 
-If `warn` is true then scaling will raise a warning.
+~~If `warn` is true then scaling will raise a warning.~~
+
+`invalid_fill_value` specifies the value to set the total cell content to if the total is less than 0
+(meaning that total tracer conservation can not be enforced). If the value is set to anything other than
+`NaN` this scheme no longer conserves mass. While this may be useful to prevent spurious numerics leading
+to crashing care should be taken that the mass doesn't deviate too much.
+
+This scheme is similar to that used by [NEMO-PISCES](https://www.nemo-ocean.eu/), although they scale the 
+tendency rather than the value, while other Earth system models simply set negative tracers to zero, for 
+example [NCAR's MARBL](https://marbl-ecosys.github.io/versions/latest_release/index.html) and
+[NEMO-TOPAZ2](https://zenodo.org/records/2648099), which does not conserve mass. More complicated schemes 
+exist, for example [ROMS-BECS](https://zenodo.org/records/3988618) uses an implicite-itterative 
+approach where each component is updated in sequence to garantee mass conservation, possibly at the 
+expense of numerical precision.
 """
-function ScaleNegativeTracers(tracers; scalefactors = ones(length(tracers)), warn = false)
+function ScaleNegativeTracers(tracers; scalefactors = ones(length(tracers)), invalid_fill_value = NaN, warn = false)
     if length(scalefactors) != length(tracers)
         error("Incorrect number of scale factors provided")
     end
 
-    return ScaleNegativeTracers(tracers, scalefactors, warn)
+    return ScaleNegativeTracers(tracers, scalefactors, invalid_fill_value, warn)
 end
 
 """
@@ -77,11 +92,11 @@ Construct a modifier to scale the conserved tracers in `model`.
 
 If `warn` is true then scaling will raise a warning.
 """
-function ScaleNegativeTracers(bgc::AbstractBiogeochemistry, grid; warn = false)
+function ScaleNegativeTracers(bgc::AbstractBiogeochemistry, grid; invalid_fill_value = NaN, warn = false)
     tracers = conserved_tracers(bgc)
     scalefactors = on_architecture(architecture(grid), ones(length(tracers)))
 
-    return ScaleNegativeTracers(tracers, scalefactors, warn)
+    return ScaleNegativeTracers(tracers, scalefactors, invalid_fill_value, warn)
 end
 
 summary(scaler::ScaleNegativeTracers) = string("Mass conserving negative scaling of $(scaler.tracers)")
@@ -97,10 +112,10 @@ function update_biogeochemical_state!(model, scale::ScaleNegativeTracers)
 
     tracers_to_scale = Tuple(model.tracers[tracer_name] for tracer_name in keys(scale.tracers))
 
-    scale_for_negs_kernel!(tracers_to_scale, scale.scalefactors)
+    scale_for_negs_kernel!(tracers_to_scale, scale.scalefactors, scale.invalid_fill_value)
 end
 
-@kernel function scale_for_negs!(tracers, scalefactors)
+@kernel function scale_for_negs!(tracers, scalefactors, invalid_fill_value)
     i, j, k = @index(Global, NTuple)
 
     t, p = 0.0, 0.0
@@ -115,7 +130,7 @@ end
         end
     end 
 
-    t < 0 && (t = NaN)
+    t < 0 && (t = invalid_fill_value)
 
     for tracer in tracers
         value = @inbounds tracer[i, j, k]