0.18.0 - Major internal refactor of NutrientsPlanktonDetritus abstraction + adds ImplicitProductivity plankton

[jagoosw]

This PR has turned into a refactor of the NPD module which hopefully aids abstraction. It also adds an ImplicitProductivity plankton (details below).

I will update docs (probably the "Implementing new models page") but this refactor hopefully aids implementing new elements by specifying the basic interfaces between nutrients, plankton, detritus, inorganic_carbon, and oxygen, but more complex models can still write their own methods up to (bgc::NutrientPlanktonDetritus{FT, <:Nutrients{<:Any, <:MyIronCycle}, <:MyNewPlankton)(i, j, k...).

This means that if you want to implement a new thing, for example a plankton, you just define the nutrient usage, elemental ratios, and waste (inorganic, dissolved, solid), and if the new model is not too abnormal it should just work. If you e.g. wanted to have a phytoplankton with variable C:N you would then need to write methods for inorganic_nitrogen_waste and inorganic_carbon_waste etc. but it would still work if you e.g. included phosphate growth limitation with no additional work.

The implicit plankton which tracks no tracers but computes a community productivity like this model from MITGCM and ClimaOceanBiogeochemistry.

This might be useful as a starting point for global runs so we can see if this simple model replicates other results (vs e.g. not being able to know if biases are from LOBTER or from the circulation, since there are no global LOBSTER solutions).

To do:

• check the examples work
• doc strings
• write docs page about the abstraction
• docs page about implicit productivity

Breaking changes:

• Moves grid and surface_PAR to be positional arguments of light models
• Changes carbonate_system to inorganic_carbon
• VariableRedfieldDetritus is now called CarbonNitrogenDissolvedParticulate

[jagoosw]

@johnryantaylor this PR maintains almost all existing functionality, except variable C:N detritus. I'm still thinking about the best way to abstract this and as its our only thing that has variable elemental ratios I think its worth considering how best to do it for e.g. when we implement MARBL plankton.

I also think that once we have devised this interface, it would be straight forward to have one of the "elements" of the large particles be calcite (or multiple calcium carbonate species/different precipitation methods), because calcite_dissolution is already a property of the detritus (and is otherwise where particulate detritus goes).

[jagoosw]

@johnryantaylor @nanophyto @lgloege

Something I've been thinking about is where to land on the balance of code complexity vs writing out similar things for different situations. As an example I implemented DissolvedParticulate to allow any number of dissolved and/or particulate classes for e.g. size or lability, which obfuscates the code to some extent because the methods gets wrapped up in expressions:

OceanBioME.jl/src/Models/AdvectedPopulations/NutrientsPlanktonDetritus/Detritus/single_element.jl

Lines 92 to 178 in 13cdab6

	function manifest_multi_class_dissolved_particulate(dissolved_names, particulate_names)
	dissolved_names = possibly_tuple_or_symbol(dissolved_names)
	particulate_names = possibly_tuple_or_symbol(particulate_names)
	for (n, name) in enumerate(dissolved_names)
	@eval begin
	@inline (bgc::NPD_DP)(i, j, k, grid, val_name::Val{$(QuoteNode(name))}, clock, fields, auxiliary_fields) = @inbounds (
	dissolved_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.dissolved_waste_partitioning[$n]
	+ dissolved_remineralisation(bgc.detritus, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.dissolved_waste_partitioning[$n]
	- grazing(bgc.plankton, bgc, i, j, k, val_name, fields, auxiliary_fields)
	- bgc.detritus.dissolved_remineralisation_rate[$n] * fields[$(QuoteNode(name))][i, j, k]
	)
	end
	end
	for (m, name) in enumerate(particulate_names)
	@eval begin
	@inline (bgc::NPD_DP)(i, j, k, grid, val_name::Val{$(QuoteNode(name))}, clock, fields, auxiliary_fields) = @inbounds (
	solid_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.particulate_waste_partitioning[$m]
	- grazing(bgc.plankton, bgc, i, j, k, val_name, fields, auxiliary_fields)
	- bgc.detritus.particulate_remineralisation_rate[$m] * fields[$(QuoteNode(name))][i, j, k]
	)
	@inline biogeochemical_drift_velocity(bgc::NPD_DP, ::Val{$(QuoteNode(name))}) =
	bgc.detritus.sinking_velocities[$m]
	end
	end
	N = length(dissolved_names)
	M = length(particulate_names)
	ex_dissolved_remineralisation = quote
	total = zero(FT)
	end
	ex_inorganic_remineralisation = quote
	total = zero(FT)
	end
	ex_calcite_remineralisation = quote
	total = zero(FT)
	end
	for (m, name) in enumerate(particulate_names)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m] *
	detritus.dissolved_fraction_of_remineralisation[$m])
	push!(ex_dissolved_remineralisation.args, ex)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m] *
	(one(FT) - detritus.dissolved_fraction_of_remineralisation[$m]))
	push!(ex_inorganic_remineralisation.args, ex)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m])
	push!(ex_calcite_remineralisation.args, ex)
	end
	for (n, name) in enumerate(dissolved_names)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.dissolved_remineralisation_rate[$n])
	push!(ex_inorganic_remineralisation.args, ex)
	end
	@eval begin
	@inline function dissolved_remineralisation(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_dissolved_remineralisation
	return total
	end
	@inline function inorganic_waste(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_inorganic_remineralisation
	return total
	end
	@inline function calcite_dissolution(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_calcite_remineralisation
	total += dissolved_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields)
	return total * carbon_ratio(bgc.plankton, bgc, i, j, k, fields) * calcite_rain_ratio(bgc.plankton, bgc, i, j, k, fields)
	end
	end
	end

I think in this case it is probably still clear enough to land on the right side of the complexity/obfuscation balance, but I've been trying to work out how to extend this for variable elemental ratio detritus and I'm not sure if it becomes too difficult to understand the code?

It gets relatively confusing because if we explicitly have the default plankton element represented (e.g. nitrogen or phosphate) only, then we have to assume that the particulates contain carbon and calcite in the same ratio as the plankton (and that the zooplankton are the same as the phytoplankton), but if we explicitly represent carbon then we can associate that calcite with the POC groups so the zooplankton doesn't have to have calcite associated with it (in LOBSTER grazing leads to some dissolution and some PIC production accounting for all the calcite in the phytoplankton). But we could also have explicit calcite groups which means the POC group now only has the organic carbon. This can be even more complicated if we then want multiple DIC/Alk because we have to have different PIC groups then because they'll dissolve differently.

(Having different DIC/Alk replicates could actually mean you need to have replicates of the plankton in some models, for example PISCES has the rain ratio in the phytoplankton vary depending on saturation state. MARBL doesn't do that and only varies rain ratio by other environmental conditions. Perhaps its not a first order effect but in the context of OAE it might actually be something we would want to capture, but would make the modelling significantly harder.)

I was wondering what you all thought of where on the code complexity/repetition we should land and if we're going to invest in more complex code that automates model setup how we should allow for future development of new elements e.g. aggregation? Perhaps there is a better abstraction for detritus that could solve some of these problems.

@glwagner do you have any thoughts on this too?

[glwagner]

@johnryantaylor @nanophyto @lgloege

Something I've been thinking about is where to land on the balance of code complexity vs writing out similar things for different situations. As an example I implemented DissolvedParticulate to allow any number of dissolved and/or particulate classes for e.g. size or lability, which obfuscates the code to some extent because the methods gets wrapped up in expressions:

OceanBioME.jl/src/Models/AdvectedPopulations/NutrientsPlanktonDetritus/Detritus/single_element.jl

Lines 92 to 178 in 13cdab6

	function manifest_multi_class_dissolved_particulate(dissolved_names, particulate_names)
	dissolved_names = possibly_tuple_or_symbol(dissolved_names)
	particulate_names = possibly_tuple_or_symbol(particulate_names)
	for (n, name) in enumerate(dissolved_names)
	@eval begin
	@inline (bgc::NPD_DP)(i, j, k, grid, val_name::Val{$(QuoteNode(name))}, clock, fields, auxiliary_fields) = @inbounds (
	dissolved_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.dissolved_waste_partitioning[$n]
	+ dissolved_remineralisation(bgc.detritus, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.dissolved_waste_partitioning[$n]
	- grazing(bgc.plankton, bgc, i, j, k, val_name, fields, auxiliary_fields)
	- bgc.detritus.dissolved_remineralisation_rate[$n] * fields[$(QuoteNode(name))][i, j, k]
	)
	end
	end
	for (m, name) in enumerate(particulate_names)
	@eval begin
	@inline (bgc::NPD_DP)(i, j, k, grid, val_name::Val{$(QuoteNode(name))}, clock, fields, auxiliary_fields) = @inbounds (
	solid_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields) * bgc.detritus.particulate_waste_partitioning[$m]
	- grazing(bgc.plankton, bgc, i, j, k, val_name, fields, auxiliary_fields)
	- bgc.detritus.particulate_remineralisation_rate[$m] * fields[$(QuoteNode(name))][i, j, k]
	)
	@inline biogeochemical_drift_velocity(bgc::NPD_DP, ::Val{$(QuoteNode(name))}) =
	bgc.detritus.sinking_velocities[$m]
	end
	end
	N = length(dissolved_names)
	M = length(particulate_names)
	ex_dissolved_remineralisation = quote
	total = zero(FT)
	end
	ex_inorganic_remineralisation = quote
	total = zero(FT)
	end
	ex_calcite_remineralisation = quote
	total = zero(FT)
	end
	for (m, name) in enumerate(particulate_names)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m] *
	detritus.dissolved_fraction_of_remineralisation[$m])
	push!(ex_dissolved_remineralisation.args, ex)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m] *
	(one(FT) - detritus.dissolved_fraction_of_remineralisation[$m]))
	push!(ex_inorganic_remineralisation.args, ex)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.particulate_remineralisation_rate[$m])
	push!(ex_calcite_remineralisation.args, ex)
	end
	for (n, name) in enumerate(dissolved_names)
	ex = :(total += fields[$(QuoteNode(name))][i, j, k] *
	detritus.dissolved_remineralisation_rate[$n])
	push!(ex_inorganic_remineralisation.args, ex)
	end
	@eval begin
	@inline function dissolved_remineralisation(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_dissolved_remineralisation
	return total
	end
	@inline function inorganic_waste(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_inorganic_remineralisation
	return total
	end
	@inline function calcite_dissolution(detritus::DissolvedParticulate{$N, $M}, bgc::NPD_DP{FT}, i, j, k, fields, auxiliary_fields) where FT
	$ex_calcite_remineralisation
	total += dissolved_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields)
	return total * carbon_ratio(bgc.plankton, bgc, i, j, k, fields) * calcite_rain_ratio(bgc.plankton, bgc, i, j, k, fields)
	end
	end
	end

I think in this case it is probably still clear enough to land on the right side of the complexity/obfuscation balance, but I've been trying to work out how to extend this for variable elemental ratio detritus and I'm not sure if it becomes too difficult to understand the code?

It gets relatively confusing because if we explicitly have the default plankton element represented (e.g. nitrogen or phosphate) only, then we have to assume that the particulates contain carbon and calcite in the same ratio as the plankton (and that the zooplankton are the same as the phytoplankton), but if we explicitly represent carbon then we can associate that calcite with the POC groups so the zooplankton doesn't have to have calcite associated with it (in LOBSTER grazing leads to some dissolution and some PIC production accounting for all the calcite in the phytoplankton). But we could also have explicit calcite groups which means the POC group now only has the organic carbon. This can be even more complicated if we then want multiple DIC/Alk because we have to have different PIC groups then because they'll dissolve differently.

(Having different DIC/Alk replicates could actually mean you need to have replicates of the plankton in some models, for example PISCES has the rain ratio in the phytoplankton vary depending on saturation state. MARBL doesn't do that and only varies rain ratio by other environmental conditions. Perhaps its not a first order effect but in the context of OAE it might actually be something we would want to capture, but would make the modelling significantly harder.)

I was wondering what you all thought of where on the code complexity/repetition we should land and if we're going to invest in more complex code that automates model setup how we should allow for future development of new elements e.g. aggregation? Perhaps there is a better abstraction for detritus that could solve some of these problems.

@glwagner do you have any thoughts on this too?

I think in general you have to take an empirical approach here... sometimes it can be better to error on the side of clarity (ie duplicating some code) at first, and then eliminate duplications if maintenance becomes hard, later.

Would it help if you develop an interface for the "tracer name type" in Oceananigans? It seems like at least some of the difficulty comes from the fact that you need to dispatch on ::Val{name}. Not sure that would help.

I noticed that you have functions like

dissolved_waste(bgc.plankton, bgc, i, j, k, fields, auxiliary_fields)

I recommend always using the "kernel form" dissolved_waste(i, j, k, grid, args...). Otherwise, users cannot write diagnostics (without a wrapper) which evaluates that RHS term. If you use the kernel form, you can call the functions from KernelFunctionOperation.

[johnryantaylor]

Should there be a - sign in front of this line? uptake of nitrate should increase the alkalinity. In the nitrate_primary_production function, nutrient_uptake is negative. Hence, without a minus sign here, it looks like nitrate uptake is reducing alkalinity instead of increasing it.

[jagoosw]

I think this is correct: #365 (comment)

[johnryantaylor]

what about the minus sign here on line 17 of defaults.jl:
- nutrient_uptake(i, j, k, grid, bgc.plankton, bgc, fields, auxiliary_fields)
I was concerned that nitrate uptake is decreasing alkalinity.

[jagoosw]

Oh I understand now, I think I've switched the nitrate and ammonia places

[johnryantaylor]

Hi @jagoosw, what are we still missing for this PR? Do you want me to work on the docs or an example?

[jagoosw]

Hi @jagoosw, what are we still missing for this PR? Do you want me to work on the docs or an example?

I'll go through and fix the docs etc in the next few days