diff --git a/Project.toml b/Project.toml index 7ad17f7c..682bbd52 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ CUDA = "5" Enzyme = "0.13, 0.14" ForwardDiff = "1" JSON = "1" -OceanBioME = "0.16, 0.18" +OceanBioME = "0.19" Oceananigans = "0.101.1, 0.102, 0.105, 0.106, 0.107, 0.108, 0.109, 0.110" SciMLBase = "2" julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index 6e1f3ee9..f54271e7 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,7 +18,7 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Documenter = "~1.17.0" ForwardDiff = "1" OrdinaryDiffEq = "6" -OceanBioME = "0.16, 0.18" +OceanBioME = "0.19" Oceananigans = "0.101.1, 0.102, 0.105, 0.106, 0.107, 0.108, 0.109, 0.110" julia = "1.10" diff --git a/docs/src/api.md b/docs/src/api.md index 379ba29f..96546ae7 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -93,13 +93,16 @@ identify a registered family and its version rather than serializing process or Agate.Processes.AbstractFormulation Agate.Processes.Smith Agate.Processes.Geider +Agate.Processes.ExponentialSaturation Agate.Processes.Monod +Agate.Processes.InhibitedMonod Agate.Processes.NormalizedDroop Agate.Processes.QuotaRegulatedMonod Agate.Processes.Liebig Agate.Processes.FrankTNorm Agate.Processes.Q10 Agate.Processes.PreferentialGrazing +Agate.Processes.LinearGrazing Agate.Processes.HeterotrophicConsumption Agate.Processes.LinearMortality Agate.Processes.QuadraticMortality @@ -136,7 +139,8 @@ The keyed parameter block separates runtime process parameters from construction Scientific slots and realized process applicability determine `Parameter` vector or matrix storage automatically, so runtime parameters never restate axes. `ConstructionParameter` values exist only during construction to feed `DerivedDefault` calculations; shaped construction parameters use the global -`axes=:plankton` construction domain. Scientific slot-to-parameter relationships are authored +`axes=:plankton` construction domain. Family-level scientific properties that affect the realized model but are not inputs to an +individual process equation use `ModelSetting` instead. Scientific slot-to-parameter relationships are authored beside the process or factor through `bindings=`. ```@docs @@ -148,7 +152,7 @@ Agate.Parameters.derive_default Agate.Parameters.NoDefault Agate.Parameters.DiameterIndexedVectorDefault Agate.Parameters.AllometricPalatability -Agate.Parameters.ConsumerAssimilation +Agate.Parameters.ConsumerResourceFromConsumer ``` ### Custom process and factor extension @@ -182,11 +186,13 @@ Agate.Compilation.process_parameter_operands ## Named families, recipes, and replay -Named model families add stable code identity and durable recipe replay around the same definition-driven process compiler. `ModelRecipe` is the `agate.model_recipe.v1` family/version/realization document: it records the registered family, exact `definition_version`, canonical plankton/size realization, parameter overrides, sinking choices, and bottom state. Named scientific mappings are serialized as mappings, so key insertion order does not change recipe equality or the scientific content hash. The loaded family supplies the canonical component/process definition on replay. `ModelManifest` records the resolved execution state. +Named model families add stable code identity and durable recipe replay around the same definition-driven process compiler. `ModelRecipe` is the `agate.model_recipe.v0.2` family/version/realization document: it records the registered family, exact `definition_version`, canonical plankton/size realization, process parameter overrides, model-setting overrides, sinking choices, and bottom state. Named scientific mappings are serialized as mappings, so key insertion order does not change recipe equality or the scientific content hash. The loaded family supplies the canonical component/process definition and setting defaults on replay. `ModelManifest` records the resolved execution state, including fully resolved model settings. -External family packages subtype `AbstractModelFamily`, provide `default_components`, `default_processes`, `definition_version`, and `parameter_definitions`, and register durable recipe identity through `family_id` and `registered_family`. Their user-facing constructors translate family-specific keywords into the nested `plankton_pfts` mapping and parameter overrides, then call `Construction.construct(family; ...)`. `normalize_pft_size_structure` provides the shared named-family `(n=0,)` shorthand without weakening core size validation. Recipes are captured with `capture_model_recipe`, the durable schema identifier is available through `recipe_schema`, and replay uses `construct(recipe)` or `construct_plus_manifest(recipe)`. +External family packages subtype `AbstractModelFamily`, provide `default_components`, `default_processes`, `definition_version`, and `parameter_definitions`, and register durable recipe identity through `family_id` and `registered_family`. Family-level scientific properties that are not inputs to an individual process equation are declared separately with `setting_definitions` and `ModelSetting`; examples include elemental or diagnostic ratios used by an integration layer. User-facing constructors translate family-specific keywords into `plankton_pfts`, process `parameter_overrides`, and optional `setting_overrides`, then call `Construction.construct(family; ...)`. `normalize_pft_size_structure` provides the shared named-family `(n=0,)` shorthand without weakening core size validation. Recipes are captured with `capture_model_recipe`, the durable schema identifier is available through `recipe_schema`, and replay uses `construct(recipe)` or `construct_plus_manifest(recipe)`. ```@docs +Agate.ModelFamilies.ModelSetting +Agate.ModelFamilies.setting_definitions Agate.Construction.ModelRecipe Agate.Construction.ModelManifest Agate.Construction.construct_plus_manifest @@ -195,6 +201,7 @@ Agate.Construction.recipe_schema Agate.Construction.normalize_pft_size_structure Agate.Construction.replay_family Agate.Construction.resolve_construction_scalar_type +Agate.Introspection.model_settings Agate.Construction.family_id Agate.Construction.registered_family Agate.Construction.encode_recipe diff --git a/docs/src/architecture_overview.md b/docs/src/architecture_overview.md index 0798edcb..f70559a9 100644 --- a/docs/src/architecture_overview.md +++ b/docs/src/architecture_overview.md @@ -13,7 +13,7 @@ A process-defined model moves through five stages: Agate validates scientific references and structural compatibility, canonicalizes process identity and factor order, canonicalizes intrinsic or family plankton realization once before layout construction, realizes each PFT into one or more SizeClasses separately from concrete prognostic tracers, discovers process drivers, and resolves process-local participant axes. 3. **Flux compilation** (`Compilation/`). - A setup-time compile context carries the canonical definition, realized layout, and parameter plan through lowering. Parameter operands resolve realized entity identity directly against storage labels during setup, so only final static indices enter runtime terms. Participant states and one-tracer-per-entity components are realized through one shared tracer + entity-position traversal, which is reused by one-axis and two-axis process lowering. Growth authoring declares `reference_resource` and any element-keyed `additional_resources` explicitly, so canonicalization no longer infers material transfer from nutrient-factor topology. Growth may have zero or more multiplicative factors; with none, its rate is simply the maximum-rate scale times biomass. Additional stoichiometric draws are valid only for Elements that the Growth plankton does not already carry as explicit prognostic States. `NutrientLimitation` may therefore mix external `NutrientResponse` and internal `QuotaResponse` subfactors while material transfer remains process-owned. Growth also owns its maximum-rate binding; Smith and Geider light factors consume that resolved process scale rather than declaring a second parameter slot. Each named process produces generic target + rate + weight flux specifications, which are grouped by target tracer and lowered into static compiled equations with no target symbols or process metadata in runtime terms. + A setup-time compile context carries the canonical definition, realized layout, and parameter plan through lowering. Parameter operands resolve realized entity identity directly against storage labels during setup, so only final static indices enter runtime terms. Participant states and one-tracer-per-entity components are realized through one shared tracer + entity-position traversal, which is reused by one-axis and two-axis process lowering. Growth authoring declares `reference_resource` and any element-keyed `additional_resources` explicitly, so canonicalization no longer infers material transfer from nutrient-factor topology. Growth may have zero or more multiplicative factors; with none, its rate is simply the maximum-rate scale times biomass. Additional stoichiometric draws are valid only for Elements that the Growth plankton does not already carry as explicit prognostic States. `NutrientLimitation` may therefore mix external `NutrientResponse` and internal `QuotaResponse` subfactors while material transfer remains process-owned. Growth also owns its maximum-rate binding; Smith and Geider light factors consume that resolved process scale rather than declaring a second parameter slot, while light formulations that do not depend on the growth scale use only their own inputs and parameters. Each named process produces generic target + rate + weight flux specifications, which are grouped by target tracer and lowered into static compiled equations with no target symbols or process metadata in runtime terms. 4. **Construction and replay** (`Construction/`). Direct `ModelDefinition` construction resolves defaults and overrides from the model @@ -84,7 +84,7 @@ Components describe structure rather than ecological role. A `Plankton` declares A PFT is defined by its functional parameter identity, not by size. Every PFT realizes at least one SizeClass. Without an explicit size structure it realizes one implicit singleton SizeClass named by the PFT, with no diameter metadata; an explicit size structure realizes one or more `_` SizeClasses with physical diameters. `n=0` is only a high-level named-model constructor shorthand for that implicit singleton and is normalized to `nothing` before core realization; `n=1` denotes one explicit SizeClass with a real diameter. -A mixotroph is an ordinary plankton participating in both growth and living-prey consumption. `Consumption` is the single consumer-resource process for living prey, bacterivory, mixotrophy, and material-pool consumption. `PreferentialGrazing` interprets `maximum_rate` as one consumer-level ingestion capacity shared across all declared living prey, while `HeterotrophicConsumption` shares one consumer-level uptake capacity across substitutable substrates according to their half-saturation and `substrate_preference`. Process products are expressed directly through `products=` or `unassimilated_products=`. `Products` provides conservative named allocation when one process flux has multiple destinations. Collection-valued participant roles use plural keywords such as `plankton=`, `consumers=`, `resources=`, and `sources=`; each accepts either one `Symbol` or a tuple and is canonicalized to a tuple during authoring. Remineralization maps many `sources=` to one `destination=`. Bacterioplankton may consume POM and be consumed as living prey through the same consumer-resource machinery. Pools are scalar material inventories; size/PFT realization is owned by `Plankton`. +A mixotroph is an ordinary plankton participating in both growth and living-prey consumption. `Consumption` is the single consumer-resource process for living prey, bacterivory, mixotrophy, and material-pool consumption. `PreferentialGrazing` interprets `maximum_rate` as one consumer-level ingestion capacity shared across all declared living prey, while `HeterotrophicConsumption` shares one consumer-level uptake capacity across substitutable substrates, with consumer-resource `half_saturation` and `substrate_preference` controlling affinity and accessibility. Process products are expressed directly through `products=` or `unassimilated_products=`. `Products` provides conservative named allocation when one process flux has multiple destinations. Collection-valued participant roles use plural keywords such as `plankton=`, `consumers=`, `resources=`, and `sources=`; each accepts either one `Symbol` or a tuple and is canonicalized to a tuple during authoring. Remineralization maps many `sources=` to one `destination=`. Bacterioplankton may consume organic-matter pools and be consumed as living prey through the same consumer-resource machinery. Pools are scalar material inventories; size/PFT realization is owned by `Plankton`. Named factors are multiplicative within a process, while independent named processes add through their fluxes to a tracer equation. Growth bookkeeping is validated per Element: an Element represented by an explicit prognostic plankton State is updated by an explicit state-changing process such as `NutrientUptake` and cannot also be supplied implicitly through Growth stoichiometry, while implicit Elements may be coupled through `FixedStoichiometry`. `NutrientLimitation.responses` is keyed by Element identity, independently of whether each response reads an external Pool or an internal quota State. Built-in Growth and `NutrientUptake` do not synthesize arbitrary non-elemental prognostic states; model families that require such synthesis can currently provide it through the custom-process extension boundary. Products and stoichiometry map process rates into affected material and element pools. Multi-destination `Products` authors exactly N-1 fractions; the omitted destination receives `1 - sum(fractions)`, so routing closes conservatively without a redundant parameter. A product may target one pool directly, derive several elemental products from a one-element source through `FixedStoichiometry`, or route the actual elemental inventories of a multi-state plankton through an element-to-pool mapping. Multi-state mortality and living-prey consumption use the reference state to define the shared specific loss intensity, apply that intensity to every prognostic state, and route only states with an Element into elemental products. Non-elemental states such as chlorophyll are removed proportionally without creating a second elemental inventory. When one mortality or living-prey consumption process routes multi-element products from multiple source plankton, those sources currently must expose the same prognostic Element set; heterogeneous source Element sets are a deferred extension rather than part of the current contract. diff --git a/examples/detritus_bacteria.jl b/examples/detritus_bacteria.jl index ba5f7d29..b6cc09c8 100644 --- a/examples/detritus_bacteria.jl +++ b/examples/detritus_bacteria.jl @@ -5,7 +5,7 @@ # graze the living bacterial plankton. A Q10 factor modifies POM consumption. using Agate.Components: Plankton, Pool -using Agate.Parameters: AllometricPalatability, ConsumerAssimilation, ConstructionParameter, DerivedDefault, Parameter +using Agate.Parameters: AllometricPalatability, ConsumerResourceFromConsumer, ConstructionParameter, DerivedDefault, Parameter using Agate.Construction: construct using Agate.Introspection: auxiliary_field_names, tracer_names using Agate.Processes: @@ -40,7 +40,7 @@ processes = ( ), factors=( temperature=Temperature( - Q10(); + Q10(:consumer); bindings=( q10=:temperature_q10, reference_temperature=:reference_temperature, @@ -90,7 +90,7 @@ parameters = ( ) ), living_assimilation=Parameter( - DerivedDefault(ConsumerAssimilation(); deps=(:assimilation_efficiency,)) + DerivedDefault(ConsumerResourceFromConsumer(); deps=(:assimilation_efficiency,)) ), ) diff --git a/examples/mixotrophy.jl b/examples/mixotrophy.jl index 5d95eb5f..845edaf5 100644 --- a/examples/mixotrophy.jl +++ b/examples/mixotrophy.jl @@ -4,7 +4,7 @@ # `Plankton` that participates in both growth and grazing. using Agate.Components: Plankton, Pool -using Agate.Parameters: AllometricPalatability, ConsumerAssimilation, ConstructionParameter, DerivedDefault, Parameter +using Agate.Parameters: AllometricPalatability, ConsumerResourceFromConsumer, ConstructionParameter, DerivedDefault, Parameter using Agate.Construction: construct using Agate.Introspection: auxiliary_field_names, tracer_names using Agate.Processes: @@ -69,7 +69,7 @@ parameters = ( ) ), assimilation_matrix=Parameter( - DerivedDefault(ConsumerAssimilation(); deps=(:assimilation_efficiency,)) + DerivedDefault(ConsumerResourceFromConsumer(); deps=(:assimilation_efficiency,)) ), ) diff --git a/scripts/lobster3_parameter_only_setup.jl b/scripts/lobster3_parameter_only_setup.jl new file mode 100644 index 00000000..aca871b4 --- /dev/null +++ b/scripts/lobster3_parameter_only_setup.jl @@ -0,0 +1,83 @@ +# Parameter-only FrankenLOBSTER approximation of the supplied LOBSTER3 model. +# +# FrankenLOBSTER defaults already match the LOBSTER3 3-D allometries for: +# P maximum growth, P nitrate affinity, H maximum uptake, H DOM affinity, +# P/Z/H mortality, Z/H assimilation, C:N, and chlorophyll:N. +# Only deliberate departures from those defaults are specified below. +# +# Temperature experiments: +# :temperature_off => Q10(P1, P2) = (1.0, 1.0) +# :p1_temperature => Q10(P1, P2) = (1.88, 1.0) +# +# Remaining functional-form approximations: +# - LOBSTER3 mass-action grazing is approached with K_G >> prey and +# g_max(d) = K_G * 15.9 * V(d)^(-0.16) / day. +# - The inactive LOBSTER3 NH4 growth branch is made negligible with a very +# large NH4 half-saturation. +# - LOBSTER3 Monod light (K=55 W m^-2) is matched at its 50% point by the +# FrankenLOBSTER exponential-saturation light response. + +using Agate +using OceanBioME +using Oceananigans +using Oceananigans.Units: day +using Agate.Library.Allometry: AllometricParam, ConstantParam, PowerLaw +using OceanBioME.Models.NutrientsPlanktonDetritusModels: DissolvedParticulate, LOBSTER + +const FrankenLOBSTER = Agate.Models.FrankenLOBSTER + +const GRAZING_K = 100.0 # >99% of mass-action rate for total prey < 1 mmol N m^-3 +const NH4_OFF_K = 1.0e12 +const LIGHT_K = 55 / log(2) # exponential response = 0.5 at PAR = 55 W m^-2 + +function temperature_q10(experiment::Symbol) + experiment === :temperature_off && return (P_1=1.0, P_2=1.0) + experiment === :p1_temperature && return (P_1=1.88, P_2=1.0) + throw(ArgumentError("temperature_experiment must be :temperature_off or :p1_temperature")) +end + +const BASE_PARAMETERS = ( + # LOBSTER3 currently uses nitrate-supported P growth only. + ammonia_half_saturation=ConstantParam(NH4_OFF_K), + nitrate_ammonia_inhibition=0.0, + + # LOBSTER3 uses Monod(PAR; K=55); FrankenLOBSTER uses exponential saturation. + light_half_saturation=(P_1=LIGHT_K, P_2=LIGHT_K), + + # Disabled in the supplied LOBSTER3 setup. + phytoplankton_exudation_fraction=(P_1=0.0, P_2=0.0), + zooplankton_excretion_rate=(Z_1=0.0, Z_2=0.0), + + # Z1 -> P1 + H1; Z2 -> P2. Scale g_max with K_G so the Holling response + # approaches LOBSTER3's mass-action g(d) * prey * predator formulation. + maximum_predation_rate=AllometricParam( + PowerLaw(); prefactor=GRAZING_K * 15.9 / day, exponent=-0.16 + ), + grazing_half_saturation=(Z_1=GRAZING_K, Z_2=GRAZING_K), + palatability_matrix=[1.0 0.0 1.0; 0.0 1.0 0.0], +) + +lobster3_parameters(; temperature_experiment=:temperature_off) = merge( + BASE_PARAMETERS, + (temperature_q10=temperature_q10(temperature_experiment),), +) + +function lobster3_like_bgc(grid; temperature_experiment=:temperature_off) + plankton = FrankenLOBSTER.construct(; + grid, + parameters=lobster3_parameters(; temperature_experiment), + ) + + # The only LOBSTER detritus default changed by the supplied LOBSTER3 setup. + detritus = DissolvedParticulate( + grid; + dissolved_remineralisation_rate=0.0, + ) + + return LOBSTER(grid; plankton, detritus) +end + +# Directly runnable one-cell setups for the two intended experiments. +grid = RectilinearGrid(CPU(); size=(1, 1, 1), extent=(1, 1, 1)) +bgc_temperature_off = lobster3_like_bgc(grid; temperature_experiment=:temperature_off) +bgc_p1_temperature = lobster3_like_bgc(grid; temperature_experiment=:p1_temperature) diff --git a/src/Agate.jl b/src/Agate.jl index 0fd3d44d..aaebd7f1 100644 --- a/src/Agate.jl +++ b/src/Agate.jl @@ -11,6 +11,7 @@ include("Runtime/Runtime.jl") include("Compilation/Compilation.jl") include("Diagnostics/Diagnostics.jl") include("Construction/Construction.jl") +include("Integrations/Integrations.jl") include("Models/Models.jl") include("Introspection.jl") @@ -25,6 +26,7 @@ export Processes export Runtime export Diagnostics export Construction +export Integrations export Introspection export ModelDefinition diff --git a/src/Compilation/Compilation.jl b/src/Compilation/Compilation.jl index ed473799..58dad963 100644 --- a/src/Compilation/Compilation.jl +++ b/src/Compilation/Compilation.jl @@ -11,7 +11,7 @@ using ..Components: using ..Processes: AbstractFactor, Growth, - Light, + Light, Smith, Geider, QuotaResponse, Consumption, Mortality, @@ -26,6 +26,7 @@ using ..Processes: CanonicalProcess, process_id, PreferentialGrazing, + LinearGrazing, Products, CanonicalModelDefinition, formulation, diff --git a/src/Compilation/consumption.jl b/src/Compilation/consumption.jl index 381d096c..2b1fd810 100644 --- a/src/Compilation/consumption.jl +++ b/src/Compilation/consumption.jl @@ -27,6 +27,27 @@ function _consumption_rate( return RateOp(formulation, operands; factors=rate_factors) end +function _consumption_rate( + formulation::LinearGrazing, + slots, + context::CompileContext, + named::CanonicalProcess, + inventory::Symbol, + _reference_resource::Symbol, + consumer::Symbol, + axis_positions::NamedTuple, + _shared_operands::Tuple, +) + operands = ( + input_operand(context.layout, inventory), + input_operand(context.layout, consumer), + parameter_operand(slots.rate, context, axis_positions), + parameter_operand(slots.palatability, context, axis_positions), + ) + rate_factors = _factor_ops(context, named, axis_positions) + return RateOp(formulation, operands; factors=rate_factors) +end + function _consumption_rate( formulation::HeterotrophicConsumption, slots, @@ -218,6 +239,13 @@ function process_fluxes( ) end end + elseif form isa LinearGrazing + for consumer in consumers, resource in resources + axis_positions = (consumer=consumer.position, resource=resource.position) + _living_consumption_fluxes!( + fluxes, named, context, consumer, resource, slots, axis_positions, () + ) + end else for consumer in consumers resource_operands = Tuple( diff --git a/src/Compilation/factors.jl b/src/Compilation/factors.jl index a861644a..d0c00aa2 100644 --- a/src/Compilation/factors.jl +++ b/src/Compilation/factors.jl @@ -35,11 +35,11 @@ function _factor_inputs(factor::QuotaResponse, named::CanonicalProcess) end function _factor_process_operands( - ::Light, + ::Light{Formulation}, context::CompileContext, named::CanonicalProcess, axis_positions::NamedTuple, -) +) where {Formulation<:Union{Smith,Geider}} ref = named.binding_refs.process.maximum_rate return (parameter_operand(ref, context, axis_positions),) end diff --git a/src/Compilation/growth.jl b/src/Compilation/growth.jl index 9d7b27e9..04788763 100644 --- a/src/Compilation/growth.jl +++ b/src/Compilation/growth.jl @@ -52,7 +52,26 @@ function process_fluxes( for participant in participants rate = _growth_rate(named, context, participant, scale_ref) - push!(fluxes, FluxSpec(participant.tracer, rate, Weight{1}())) + product_targets = named.semantic_facts.product_targets + if isnothing(product_targets) + push!(fluxes, FluxSpec(participant.tracer, rate, Weight{1}())) + else + axis_positions = (plankton=participant.position,) + product_fraction = parameter_operand( + named.binding_refs.process.product_fraction, context, axis_positions + ) + retained_fraction = ComplementOp((product_fraction,)) + push!( + fluxes, + FluxSpec(participant.tracer, rate, Weight{1}((retained_fraction,))), + ) + append!( + fluxes, + _product_fluxes( + named, product_targets, context, rate; suffix=(product_fraction,) + ), + ) + end append!(fluxes, _growth_resource_fluxes(named, context, rate)) end return Tuple(fluxes) diff --git a/src/Components/layout.jl b/src/Components/layout.jl index 076bb7f8..986f9b3b 100644 --- a/src/Components/layout.jl +++ b/src/Components/layout.jl @@ -373,6 +373,7 @@ function model_metadata(layout::ModelLayout; parameter_axes=(;), parameter_const ) return (; pft_entities, + component_tracers=layout.component_tracers, plankton_tracers, plankton_diameters=Tuple(diameter_metadata(diameter) for diameter in layout.size_class_diameters), parameter_axes, diff --git a/src/Construction/Construction.jl b/src/Construction/Construction.jl index 82b80fc9..83e31d2d 100644 --- a/src/Construction/Construction.jl +++ b/src/Construction/Construction.jl @@ -2,12 +2,14 @@ module Construction import Oceananigans -export construct, construct_plus_manifest +export construct, construct_plus_manifest, construct_plus_recipe export ModelRecipe, ModelManifest export capture_model_recipe, recipe_schema -export normalize_pft_size_structure +export normalize_pft_size_structure, plankton_realization +export resolve_model_settings export encode_recipe, decode_recipe, export_recipe, import_recipe +include("settings.jl") include("recipe.jl") include("recipe_serialization.jl") include("recipe_provenance.jl") diff --git a/src/Construction/construct.jl b/src/Construction/construct.jl index 2f00f600..0e9e42bf 100644 --- a/src/Construction/construct.jl +++ b/src/Construction/construct.jl @@ -6,7 +6,7 @@ import Oceananigans using Oceananigans.Architectures: architecture, CPU, GPU -using ..ModelFamilies: AbstractModelFamily +using ..ModelFamilies: AbstractModelFamily, default_components, plankton_roles using ..Components: canonicalize_plankton_realization, realize_model_layout, model_metadata @@ -112,6 +112,28 @@ end return n isa Integer && !(n isa Bool) && n == 0 ? nothing : specification end +"""Translate a family's user-facing plankton roles into canonical logical PFT realization.""" +function plankton_realization(family::AbstractModelFamily, size_structure) + size_structure isa NamedTuple || throw(ArgumentError("size_structure must be a NamedTuple")) + roles = plankton_roles(family) + role_names = keys(roles) + Set(keys(size_structure)) == Set(role_names) || throw( + ArgumentError("size_structure must define exactly $(collect(role_names))") + ) + + component_names = Tuple(values(roles)) + component_values = ntuple(length(role_names)) do i + role = role_names[i] + pfts = getproperty(size_structure, role) + pfts isa NamedTuple || throw( + ArgumentError("size_structure.$role must be a NamedTuple") + ) + NamedTuple{keys(pfts)}(Tuple(normalize_pft_size_structure(value) for value in values(pfts))) + end + authored = NamedTuple{component_names}(component_values) + return canonicalize_plankton_realization(default_components(family), authored) +end + function _realize_process_definition( definition, @@ -129,6 +151,7 @@ function _construct_process_definition( definition::ModelDefinition; plankton_pfts=nothing, parameter_overrides::NamedTuple=(;), + model_settings::NamedTuple=(;), sinking_tracers=nothing, open_bottom::Bool=true, grid=nothing, @@ -203,10 +226,13 @@ function _construct_process_definition( runtime_parameters = runtime_parameter_values(parameter_plan, resolved_parameters) compile_context = CompileContext(canonical, layout, parameter_plan) equations = compile_model_tendencies(compile_context; target_order=tracer_names) - metadata = model_metadata( - layout; - parameter_axes=parameter_plan_metadata(canonical, parameter_plan), - parameter_constraints=constraints, + metadata = merge( + model_metadata( + layout; + parameter_axes=parameter_plan_metadata(canonical, parameter_plan), + parameter_constraints=constraints, + ), + (; model_settings), ) sinking_velocities = isnothing(sinking_tracers) ? nothing : setup_velocity_fields(sinking_tracers, grid, open_bottom) @@ -221,6 +247,7 @@ function _construct_process_definition( capture_model_manifest( manifest_family, resolved_parameters, + model_settings, layout, parameter_plan; tracer_order=tracer_names, @@ -245,9 +272,18 @@ function _construct_registered_model( scalar_type=nothing, build_manifest::Bool=false, ) + T = resolve_construction_scalar_type(grid, scalar_type) + model_settings = resolve_model_settings(family, realization.setting_overrides, T) + process_realization = (; + plankton_pfts=realization.plankton_pfts, + parameter_overrides=realization.parameter_overrides, + sinking_tracers=realization.sinking_tracers, + open_bottom=realization.open_bottom, + ) return _construct_process_definition( ModelDefinition(family); - realization..., + process_realization..., + model_settings, grid, arch, scalar_type, @@ -264,9 +300,11 @@ function _construct_recipe( scalar_type=nothing, build_manifest::Bool=false, ) + family = replay_family(recipe) + realization = _family_realization(recipe) return _construct_registered_model( - replay_family(recipe), - _family_realization(recipe); + family, + realization; grid, arch, scalar_type, @@ -277,28 +315,33 @@ end """ construct(family::AbstractModelFamily; - plankton_pfts, parameter_overrides=(;), + plankton_pfts, parameter_overrides=(;), setting_overrides=(;), sinking_tracers=nothing, open_bottom=true, grid=nothing, arch=nothing, scalar_type=nothing) -> bgc Construct a registered model family from its resolved family realization. This is the supported construction seam for external family packages after their own user-facing -constructor syntax has been translated into the nested `plankton_pfts` mapping and -parameter overrides. Runtime grid, architecture, and scalar precision remain execution -choices. +constructor syntax has been translated into the nested `plankton_pfts` mapping, process +parameter overrides, and family-level setting overrides. Runtime grid, architecture, and scalar +precision remain execution choices. """ function construct( family::AbstractModelFamily; plankton_pfts::NamedTuple, parameter_overrides::NamedTuple=(;), + setting_overrides::NamedTuple=(;), sinking_tracers=nothing, open_bottom::Bool=true, grid=nothing, arch=nothing, scalar_type=nothing, ) - realization = (; plankton_pfts, parameter_overrides, sinking_tracers, open_bottom) - bgc, _ = _construct_registered_model(family, realization; grid, arch, scalar_type) + realization = (; + plankton_pfts, parameter_overrides, setting_overrides, sinking_tracers, open_bottom + ) + bgc, _ = _construct_registered_model( + family, realization; grid, arch, scalar_type + ) return bgc end @@ -342,12 +385,34 @@ end """Replay a versioned family recipe in the supplied execution environment.""" function construct( - recipe::ModelRecipe; grid=nothing, arch=nothing, scalar_type=nothing + recipe::ModelRecipe; + grid=nothing, + arch=nothing, + scalar_type=nothing, ) bgc, _ = _construct_recipe(recipe; grid, arch, scalar_type) return bgc end +"""Construct a registered family and capture the canonical recipe used for construction.""" +function construct_plus_recipe( + family::AbstractModelFamily; + plankton_pfts::NamedTuple, + parameter_overrides::NamedTuple=(;), + setting_overrides::NamedTuple=(;), + sinking_tracers=nothing, + open_bottom::Bool=true, + grid=nothing, + arch=nothing, + scalar_type=nothing, +) + recipe = capture_model_recipe( + family; plankton_pfts, parameter_overrides, setting_overrides, sinking_tracers, open_bottom + ) + bgc = construct(recipe; grid, arch, scalar_type) + return bgc, recipe +end + """Replay a versioned family recipe and return its resolved manifest.""" function construct_plus_manifest( recipe::ModelRecipe; grid=nothing, arch=nothing, scalar_type=nothing diff --git a/src/Construction/parameter_realization.jl b/src/Construction/parameter_realization.jl index 75d3d676..69462fb4 100644 --- a/src/Construction/parameter_realization.jl +++ b/src/Construction/parameter_realization.jl @@ -153,15 +153,16 @@ function materialize_parameter_default( provider::ConstantDefault, parameter, ::Type{T} ) where {T<:Real} value = provider.value - value = value isa Bool ? value : T(value) rank = parameter.rank - rank == 0 && return value + rank == 0 && return value isa Bool ? value : T(value) expected = parameter.storage_shape if rank == 1 - return fill(value, only(expected)) + value isa AbstractVector && return materialize_parameter_value(parameter, value, T) + return fill(value isa Bool ? value : T(value), only(expected)) elseif rank == 2 - return fill(value, expected...) + value isa AbstractMatrix && return materialize_parameter_value(parameter, value, T) + return fill(value isa Bool ? value : T(value), expected...) end throw(ArgumentError("parameter :$(parameter.name) has unsupported rank $rank")) end diff --git a/src/Construction/recipe.jl b/src/Construction/recipe.jl index f69cf01d..da50ac9c 100644 --- a/src/Construction/recipe.jl +++ b/src/Construction/recipe.jl @@ -15,28 +15,32 @@ end """Versioned registered-family recipe captured before runtime realization. `ModelRecipe` stores only the registered family identity, its exact scientific -`definition_version`, and canonical construction inputs. `==`, `isequal`, `hash`, and content +`definition_version`, and canonical construction inputs, including separate process-parameter +and model-setting overrides. `==`, `isequal`, `hash`, and content hashing share that scientific identity; named mapping insertion order is ignored. Components, processes, parameter definitions, runtime precision, host fields, and compiled equations are supplied by the loaded family implementation on replay. """ -struct ModelRecipe{PlanktonPFTs,ParameterOverrides,SinkingTracers} +struct ModelRecipe{PlanktonPFTs,ParameterOverrides,SettingOverrides,SinkingTracers} family::Symbol definition_version::VersionNumber plankton_pfts::PlanktonPFTs parameter_overrides::ParameterOverrides + setting_overrides::SettingOverrides sinking_tracers::SinkingTracers open_bottom::Bool end """Resolved deterministic scientific state produced by model construction. -`ModelManifest` records the fully materialized parameters, realized PFT entities and -tracer ordering, interaction sources, sinking configuration, and scalar type. Equality and hashing use this +`ModelManifest` records the fully materialized process parameters and model settings, realized +PFT entities and tracer ordering, interaction sources, sinking configuration, and scalar type. +Equality and hashing use this resolved scientific content; durable replay is defined by the corresponding recipe representation. """ struct ModelManifest{ Parameters, + Settings, PFTEntities, TracerOrder, AuxiliaryFields, @@ -46,6 +50,7 @@ struct ModelManifest{ ScalarType<:Real, } parameters::Parameters + settings::Settings pft_entities::PFTEntities tracer_order::TracerOrder auxiliary_fields::AuxiliaryFields @@ -65,6 +70,7 @@ _recipe_identity(recipe::ModelRecipe) = _recipe_identity( _manifest_identity(manifest::ModelManifest) = (; parameters=manifest.parameters, + settings=manifest.settings, pft_entities=manifest.pft_entities, tracer_order=manifest.tracer_order, auxiliary_fields=manifest.auxiliary_fields, @@ -96,6 +102,7 @@ function capture_model_recipe( family::AbstractModelFamily; plankton_pfts::NamedTuple, parameter_overrides::NamedTuple=(;), + setting_overrides::NamedTuple=(;), sinking_tracers=nothing, open_bottom::Bool=true, ) @@ -113,6 +120,7 @@ function capture_model_recipe( version, deepcopy(plankton_pfts), deepcopy(parameter_overrides), + deepcopy(setting_overrides), deepcopy(sinking_tracers), open_bottom, ) @@ -121,6 +129,7 @@ end _family_realization(recipe::ModelRecipe) = (; plankton_pfts=recipe.plankton_pfts, parameter_overrides=recipe.parameter_overrides, + setting_overrides=recipe.setting_overrides, sinking_tracers=recipe.sinking_tracers, open_bottom=recipe.open_bottom, ) @@ -152,6 +161,7 @@ replay_family(recipe::ModelRecipe) = function capture_model_manifest( family::AbstractModelFamily, parameters, + settings, layout::ModelLayout, parameter_plan; tracer_order::Tuple, @@ -187,6 +197,7 @@ function capture_model_manifest( return ModelManifest( deepcopy(parameters), + deepcopy(settings), pft_entities, tracer_order, auxiliary_fields, diff --git a/src/Construction/recipe_serialization.jl b/src/Construction/recipe_serialization.jl index 3d10829b..8af2045c 100644 --- a/src/Construction/recipe_serialization.jl +++ b/src/Construction/recipe_serialization.jl @@ -7,7 +7,7 @@ using ..Library.Allometry: allometric_relationship_identifier, allometric_relationship_from_identifier -const MODEL_RECIPE_SCHEMA = "agate.model_recipe.v1" +const MODEL_RECIPE_SCHEMA = "agate.model_recipe.v0.2" """Return the durable model-recipe schema identifier supported by this Agate version.""" recipe_schema() = MODEL_RECIPE_SCHEMA @@ -15,7 +15,7 @@ const _RECIPE_DOCUMENT_KEYS = ( "schema", "family", "definition_version", "realization", "provenance", "content_hash" ) const _REALIZATION_KEYS = ( - "plankton_pfts", "parameter_overrides", "sinking_tracers", "open_bottom" + "plankton_pfts", "parameter_overrides", "setting_overrides", "sinking_tracers", "open_bottom" ) const _SUPPORTED_SPACING = (:linear, :log) @@ -347,6 +347,7 @@ function _encode_realization(recipe::ModelRecipe) return Dict{String,Any}( "plankton_pfts" => _encode_plankton_pfts(recipe.plankton_pfts), "parameter_overrides" => _encode_parameter_overrides(recipe.parameter_overrides), + "setting_overrides" => _encode_parameter_overrides(recipe.setting_overrides), "sinking_tracers" => isnothing(recipe.sinking_tracers) ? nothing : _encode_parameter_overrides(recipe.sinking_tracers), "open_bottom" => recipe.open_bottom, @@ -361,12 +362,15 @@ function _decode_realization(x, path) parameter_overrides = _decode_parameter_overrides( realization["parameter_overrides"], "$path.parameter_overrides" ) + setting_overrides = _decode_parameter_overrides( + realization["setting_overrides"], "$path.setting_overrides" + ) sinking_tracers = isnothing(realization["sinking_tracers"]) ? nothing : _decode_parameter_overrides( realization["sinking_tracers"], "$path.sinking_tracers" ) open_bottom = _boolean(realization["open_bottom"], "$path.open_bottom") - return (; plankton_pfts, parameter_overrides, sinking_tracers, open_bottom) + return (; plankton_pfts, parameter_overrides, setting_overrides, sinking_tracers, open_bottom) end """Encode a versioned family recipe with a scientific content hash and package provenance.""" @@ -412,6 +416,7 @@ function decode_recipe(document::AbstractDict) version, plankton_pfts, realization.parameter_overrides, + realization.setting_overrides, realization.sinking_tracers, realization.open_bottom, ) diff --git a/src/Construction/settings.jl b/src/Construction/settings.jl new file mode 100644 index 00000000..65f1ccff --- /dev/null +++ b/src/Construction/settings.jl @@ -0,0 +1,36 @@ +using ..ModelFamilies: AbstractModelFamily, ModelSetting, setting_definitions +using ..Processes: parameter_domain_valid + +function _typed_setting(value, ::Type{T}) where {T<:Real} + value isa Real && !(value isa Bool) && return convert(T, value) + return value +end + +"""Resolve and validate family-level scientific settings for one construction.""" +function resolve_model_settings( + family::AbstractModelFamily, overrides::NamedTuple, ::Type{T} +) where {T<:Real} + definitions = setting_definitions(family) + unknown = Tuple(name for name in keys(overrides) if !hasproperty(definitions, name)) + isempty(unknown) || throw( + ArgumentError("unknown model setting override(s): $(join(string.(unknown), ", "))"), + ) + + names = keys(definitions) + values = ntuple(length(names)) do i + name = names[i] + definition = getproperty(definitions, name) + definition isa ModelSetting || throw( + ArgumentError("setting_definitions must contain only ModelSetting values"), + ) + value = hasproperty(overrides, name) ? getproperty(overrides, name) : definition.default + value = _typed_setting(value, T) + parameter_domain_valid(value, definition.domain) || throw( + ArgumentError( + "model setting :$name must satisfy domain :$(definition.domain); got $(repr(value))" + ), + ) + value + end + return NamedTuple{names}(values) +end diff --git a/src/Integrations/Integrations.jl b/src/Integrations/Integrations.jl new file mode 100644 index 00000000..b0cc754c --- /dev/null +++ b/src/Integrations/Integrations.jl @@ -0,0 +1,10 @@ +"""Integration layers between compiled Agate runtimes and external biogeochemistry frameworks.""" +module Integrations + +include("oceanbiome_npd.jl") + +export NPDPlankton +export construct_npd_plankton, construct_npd_plankton_plus_recipe +export npd_configuration + +end # module diff --git a/src/Integrations/oceanbiome_npd.jl b/src/Integrations/oceanbiome_npd.jl new file mode 100644 index 00000000..fda31163 --- /dev/null +++ b/src/Integrations/oceanbiome_npd.jl @@ -0,0 +1,451 @@ +using Adapt: adapt +import Adapt: adapt_structure + +using ..ModelFamilies: AbstractModelFamily, plankton_roles +import ..Construction + +import OceanBioME: chlorophyll +import Oceananigans.Biogeochemistry: + biogeochemical_drift_velocity, + required_biogeochemical_auxiliary_fields, + required_biogeochemical_tracers + +using OceanBioME.Models.NutrientsPlanktonDetritusModels: NutrientsPlanktonDetritus +import OceanBioME.Models.NutrientsPlanktonDetritusModels: + carbon_ratio, + chlorophyll_ratio, + dissolved_waste, + inorganic_waste, + nutrient_uptake, + solid_waste +import OceanBioME.Models.NutrientsPlanktonDetritusModels.DetritusModels: grazing + +"""OceanBioME NPD plankton component backed by one compiled Agate runtime. + +`NPDPlankton` is an integration boundary, not a biological model. The compiled +Agate runtime owns the ecological equations; this wrapper maps their signed +tracer tendencies onto OceanBioME's existing NPD plankton hooks. +""" +struct NPDPlankton{ + Runtime, + OwnedTracers, + OwnedTracerType, + NutrientTracers, + ExchangeTracers, + ConsumedDetritus, + Dependencies, + PhytoplanktonTracers, + Traits, +} + runtime::Runtime + traits::Traits +end + +function _component_tracers(runtime, components::Tuple) + metadata = runtime.metadata + hasproperty(metadata, :component_tracers) || throw( + ArgumentError("Agate runtime metadata does not expose component tracer identities."), + ) + return Tuple( + tracer + for component in components + for tracer in begin + hasproperty(metadata.component_tracers, component) || throw( + ArgumentError("Unknown Agate component :$component."), + ) + getproperty(metadata.component_tracers, component) + end + ) +end + +function _validate_npd_traits(traits::NamedTuple) + for name in (:carbon_ratio, :chlorophyll_ratio) + hasproperty(traits, name) || throw( + ArgumentError("NPDPlankton traits must define :$name."), + ) + value = getproperty(traits, name) + value isa Real && !(value isa Bool) && isfinite(value) && value >= 0 || throw( + ArgumentError("NPDPlankton trait :$name must be finite and nonnegative."), + ) + end + traits.carbon_ratio > 0 || throw(ArgumentError("NPDPlankton carbon_ratio must be > 0.")) + return traits +end + +""" + NPDPlankton(runtime; owned_components, phytoplankton_components=(), nutrient_tracers=(), + exchange_tracers=(solid=:solid_waste, dissolved=:dissolved_waste, + inorganic=:inorganic_waste), + consumed_detritus=(), dependencies=(), traits) + +Wrap a compiled Agate runtime as an OceanBioME `NutrientsPlanktonDetritus` plankton component. +Component names are resolved once from Agate runtime metadata; all cell-level coupling is then +statically dispatched from the resulting tracer tuples. `nutrient_tracers` names runtime tracers +exposed through OceanBioME's nutrient-uptake hooks; Agate does not impose a nutrient-name +vocabulary. A dependency that also names a runtime auxiliary driver is read from the surrounding +NPD tracer fields rather than requested as a separate auxiliary field. +""" +function NPDPlankton( + runtime; + owned_components::Tuple, + phytoplankton_components::Tuple=(), + nutrient_tracers::Tuple=(), + exchange_tracers::NamedTuple=( + solid=:solid_waste, + dissolved=:dissolved_waste, + inorganic=:inorganic_waste, + ), + consumed_detritus::Tuple=(), + dependencies::Tuple=(), + traits::NamedTuple, +) + keys(exchange_tracers) == (:solid, :dissolved, :inorganic) || throw( + ArgumentError("exchange_tracers must define (:solid, :dissolved, :inorganic)."), + ) + owned = _component_tracers(runtime, owned_components) + isempty(owned) && throw(ArgumentError("NPDPlankton must own at least one tracer.")) + phytoplankton = _component_tracers(runtime, phytoplankton_components) + isempty(phytoplankton) && throw( + ArgumentError("NPDPlankton must identify at least one phytoplankton tracer."), + ) + owned_type = mapreduce(name -> typeof(Val(name)), (A, B) -> Union{A,B}, owned) + runtime_tracers = required_biogeochemical_tracers(runtime) + for tracer in nutrient_tracers + tracer isa Symbol || throw( + ArgumentError("NPDPlankton nutrient_tracers must contain tracer names as Symbols."), + ) + tracer in runtime_tracers || throw( + ArgumentError("NPDPlankton nutrient tracer :$tracer is not an Agate runtime tracer."), + ) + end + for (channel, tracer) in pairs(exchange_tracers) + tracer === nothing && continue + tracer in runtime_tracers || throw(ArgumentError( + "NPDPlankton exchange channel :$channel names unknown runtime tracer :$tracer.", + )) + end + exchanges = Tuple(values(exchange_tracers)) + _validate_npd_traits(traits) + + return NPDPlankton{ + typeof(runtime), + owned, + owned_type, + nutrient_tracers, + exchanges, + consumed_detritus, + dependencies, + phytoplankton, + typeof(traits), + }(runtime, traits) +end + +"""Return family-specific `NPDPlankton` keyword overrides for a realized runtime. + +Ownership, phytoplankton identity, external dependencies, standard exchange channels, and +resolved model settings are inferred from Agate metadata. External families only need to +declare choices that cannot be inferred safely, such as nutrient uptake tracers, consumed +detritus, or physical-driver dependencies. Configured `dependencies` are appended to the +inferred tracer dependencies. +""" +npd_configuration(::AbstractModelFamily, _runtime) = (;) + +function _standard_exchange_tracers(runtime) + tracers = required_biogeochemical_tracers(runtime) + present(name) = name in tracers ? name : nothing + return ( + solid=present(:solid_waste), + dissolved=present(:dissolved_waste), + inorganic=present(:inorganic_waste), + ) +end + +function _npd_default_configuration(family::AbstractModelFamily, runtime) + roles = plankton_roles(family) + hasproperty(roles, :phytoplankton) || throw( + ArgumentError("NPD plankton families must define a :phytoplankton role."), + ) + return (; + owned_components=Tuple(unique(values(roles))), + phytoplankton_components=(roles.phytoplankton,), + nutrient_tracers=(), + exchange_tracers=_standard_exchange_tracers(runtime), + consumed_detritus=(), + traits=runtime.metadata.model_settings, + ) +end + +function _npd_dependencies(runtime, configuration) + owned = _component_tracers(runtime, configuration.owned_components) + exchanges = Tuple( + value for value in values(configuration.exchange_tracers) if value !== nothing + ) + return Tuple( + tracer for tracer in required_biogeochemical_tracers(runtime) + if !(tracer in owned) && !(tracer in exchanges) + ) +end + +function _wrap_npd_runtime(family::AbstractModelFamily, runtime) + overrides = npd_configuration(family, runtime) + configuration = merge(_npd_default_configuration(family, runtime), overrides) + inferred = _npd_dependencies(runtime, configuration) + configured = hasproperty(overrides, :dependencies) ? overrides.dependencies : () + configuration = merge( + configuration, (; dependencies=_append_unique(inferred, configured)) + ) + return NPDPlankton(runtime; configuration...) +end + +_require_npd_grid(sinking_tracers, grid) = + !isnothing(sinking_tracers) && isnothing(grid) ? + throw(ArgumentError("grid is required when `sinking_tracers` are configured")) : nothing + +"""Construct a registered Agate family directly as an OceanBioME NPD plankton component.""" +function construct_npd_plankton( + family::AbstractModelFamily; + plankton_pfts::NamedTuple, + parameter_overrides::NamedTuple=(;), + setting_overrides::NamedTuple=(;), + sinking_tracers=nothing, + open_bottom::Bool=true, + grid=nothing, + arch=nothing, + scalar_type=nothing, +) + _require_npd_grid(sinking_tracers, grid) + runtime = Construction.construct( + family; + plankton_pfts, + parameter_overrides, + setting_overrides, + sinking_tracers, + open_bottom, + grid, + arch, + scalar_type, + ) + return _wrap_npd_runtime(family, runtime) +end + +"""Construct an NPD plankton component and capture the canonical Agate family recipe.""" +function construct_npd_plankton_plus_recipe( + family::AbstractModelFamily; + plankton_pfts::NamedTuple, + parameter_overrides::NamedTuple=(;), + setting_overrides::NamedTuple=(;), + sinking_tracers=nothing, + open_bottom::Bool=true, + grid=nothing, + arch=nothing, + scalar_type=nothing, +) + _require_npd_grid(sinking_tracers, grid) + runtime, recipe = Construction.construct_plus_recipe( + family; + plankton_pfts, + parameter_overrides, + setting_overrides, + sinking_tracers, + open_bottom, + grid, + arch, + scalar_type, + ) + return _wrap_npd_runtime(family, runtime), recipe +end + +"""Replay a registered Agate family recipe directly as an OceanBioME NPD plankton component.""" +function construct_npd_plankton( + recipe::Construction.ModelRecipe; grid=nothing, arch=nothing, scalar_type=nothing +) + family = Construction.replay_family(recipe) + _require_npd_grid(recipe.sinking_tracers, grid) + runtime = Construction.construct( + recipe; + grid, + arch, + scalar_type, + ) + return _wrap_npd_runtime(family, runtime) +end + +@inline _owned_tracers(::NPDPlankton{<:Any,OwnedTracers}) where {OwnedTracers} = OwnedTracers +@inline _nutrient_tracers(::NPDPlankton{<:Any,<:Any,<:Any,NutrientTracers}) where {NutrientTracers} = NutrientTracers +@inline _exchange_tracers(::NPDPlankton{<:Any,<:Any,<:Any,<:Any,ExchangeTracers}) where {ExchangeTracers} = ExchangeTracers +@inline _consumed_detritus(::NPDPlankton{<:Any,<:Any,<:Any,<:Any,<:Any,ConsumedDetritus}) where {ConsumedDetritus} = ConsumedDetritus +@inline _dependencies(::NPDPlankton{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,Dependencies}) where {Dependencies} = Dependencies +@inline phytoplankton_tracers(::NPDPlankton{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,PhytoplanktonTracers}) where {PhytoplanktonTracers} = PhytoplanktonTracers + +@inline required_biogeochemical_tracers(plankton::NPDPlankton) = _owned_tracers(plankton) +@inline required_biogeochemical_auxiliary_fields(plankton::NPDPlankton) = Tuple( + name for name in required_biogeochemical_auxiliary_fields(plankton.runtime) + if !(name in _dependencies(plankton)) +) +@inline biogeochemical_drift_velocity(plankton::NPDPlankton, tracer::Val) = + biogeochemical_drift_velocity(plankton.runtime, tracer) + +@inline chlorophyll_ratio(plankton::NPDPlankton) = plankton.traits.chlorophyll_ratio +@inline carbon_ratio( + plankton::NPDPlankton, ::NutrientsPlanktonDetritus{FloatType} +) where FloatType = convert(FloatType, plankton.traits.carbon_ratio) +@inline chlorophyll(plankton::NPDPlankton, model) = plankton.traits.chlorophyll_ratio * + mapreduce(name -> getproperty(model.tracers, name), +, phytoplankton_tracers(plankton)) + +@inline function adapt_structure( + to, + plankton::NPDPlankton{ + <:Any,OwnedTracers,OwnedTracerType,NutrientTracers,ExchangeTracers, + ConsumedDetritus,Dependencies,PhytoplanktonTracers,<:Any, + }, +) where { + OwnedTracers,OwnedTracerType,NutrientTracers,ExchangeTracers, + ConsumedDetritus,Dependencies,PhytoplanktonTracers, +} + runtime = adapt(to, plankton.runtime) + traits = adapt(to, plankton.traits) + return NPDPlankton{ + typeof(runtime),OwnedTracers,OwnedTracerType,NutrientTracers,ExchangeTracers, + ConsumedDetritus,Dependencies,PhytoplanktonTracers,typeof(traits), + }(runtime, traits) +end + +@inline function _append_unique(acc::Tuple, values::Tuple) + isempty(values) && return acc + head = first(values) + next = head in acc ? acc : (acc..., head) + return _append_unique(next, Base.tail(values)) +end + +@inline function required_biogeochemical_tracers( + npd::NutrientsPlanktonDetritus{<:Any,<:Any,PlanktonType}, +) where {PlanktonType<:NPDPlankton} + tracers = ( + required_biogeochemical_tracers(npd.nutrients)..., + required_biogeochemical_tracers(npd.plankton)..., + required_biogeochemical_tracers(npd.detritus)..., + required_biogeochemical_tracers(npd.inorganic_carbon)..., + required_biogeochemical_tracers(npd.oxygen)..., + _dependencies(npd.plankton)..., + ) + return _append_unique((), tracers) +end + +# OceanBioME fields -> Agate's statically ordered positional state. +@inline function _runtime_tracer_value(::Val{Tracer}, plankton::NPDPlankton, i, j, k, fields) where Tracer + Tracer in _exchange_tracers(plankton) && + return zero(@inbounds getproperty(fields, first(_owned_tracers(plankton)))[i, j, k]) + return @inbounds getproperty(fields, Tracer)[i, j, k] +end + +@inline function _runtime_tracer_values(plankton::NPDPlankton, i, j, k, fields) + tracers = required_biogeochemical_tracers(plankton.runtime) + return ntuple(Val(length(tracers))) do n + _runtime_tracer_value(Val(tracers[n]), plankton, i, j, k, fields) + end +end + +@inline function _runtime_auxiliary_value( + ::Val{Auxiliary}, plankton::NPDPlankton, i, j, k, fields, auxiliary_fields +) where Auxiliary + Auxiliary in _dependencies(plankton) && + return @inbounds getproperty(fields, Auxiliary)[i, j, k] + return @inbounds getproperty(auxiliary_fields, Auxiliary)[i, j, k] +end + +@inline function _runtime_auxiliary_values( + plankton::NPDPlankton, i, j, k, fields, auxiliary_fields +) + auxiliaries = required_biogeochemical_auxiliary_fields(plankton.runtime) + return ntuple(Val(length(auxiliaries))) do n + _runtime_auxiliary_value( + Val(auxiliaries[n]), plankton, i, j, k, fields, auxiliary_fields + ) + end +end + +@inline function _agate_tendency( + plankton::NPDPlankton, tracer::Val, i, j, k, t, fields, auxiliary_fields +) + tracer_values = _runtime_tracer_values(plankton, i, j, k, fields) + auxiliary_values = _runtime_auxiliary_values( + plankton, i, j, k, fields, auxiliary_fields + ) + x = zero(t) + return plankton.runtime(tracer, x, x, x, t, tracer_values..., auxiliary_values...) +end + +@inline _exchange_tendency(plankton, tracer, i, j, k, grid, fields, auxiliary_fields) = + _agate_tendency(plankton, tracer, i, j, k, zero(eltype(grid)), fields, auxiliary_fields) + +# Restrict the NPD call overload to the Agate-owned living tracer union so OceanBioME +# nutrient/detritus/carbon/oxygen tracers keep their native dispatch. +@inline (bgc::NutrientsPlanktonDetritus{<:Any,<:Any,PlanktonType})( + i, j, k, grid, tracer::OwnedTracerType, clock, fields, auxiliary_fields +) where { + OwnedTracerType, + PlanktonType<:NPDPlankton{<:Any,<:Any,OwnedTracerType}, +} = _agate_tendency( + bgc.plankton, tracer, i, j, k, clock.time, fields, auxiliary_fields +) + +@inline function nutrient_uptake( + i, j, k, grid, nutrient::Val{Nutrient}, plankton::NPDPlankton, + ::NutrientsPlanktonDetritus, fields, auxiliary_fields, +) where Nutrient + Nutrient in _nutrient_tracers(plankton) || return zero(eltype(grid)) + return -_exchange_tendency(plankton, nutrient, i, j, k, grid, fields, auxiliary_fields) +end + +@inline _sum_nutrient_uptake( + ::Tuple{}, i, j, k, grid, plankton, bgc, fields, auxiliary_fields, +) = zero(eltype(grid)) + +@inline function _sum_nutrient_uptake( + nutrients::Tuple, i, j, k, grid, plankton, bgc, fields, auxiliary_fields, +) + nutrient = first(nutrients) + return nutrient_uptake( + i, j, k, grid, Val(nutrient), plankton, bgc, fields, auxiliary_fields + ) + _sum_nutrient_uptake( + Base.tail(nutrients), i, j, k, grid, plankton, bgc, fields, auxiliary_fields + ) +end + +@inline function nutrient_uptake( + i, j, k, grid, plankton::NPDPlankton, + bgc::NutrientsPlanktonDetritus, fields, auxiliary_fields, +) + return _sum_nutrient_uptake( + _nutrient_tracers(plankton), i, j, k, grid, plankton, bgc, fields, auxiliary_fields + ) +end + +@inline function _exchange_channel( + plankton, channel, i, j, k, grid, fields, auxiliary_fields, +) + channel === nothing && return zero(eltype(grid)) + return _exchange_tendency( + plankton, Val(channel), i, j, k, grid, fields, auxiliary_fields + ) +end + +for (hook, index) in ((:solid_waste, 1), (:dissolved_waste, 2), (:inorganic_waste, 3)) + @eval @inline function $hook( + i, j, k, grid, plankton::NPDPlankton, + ::NutrientsPlanktonDetritus, fields, auxiliary_fields, + ) + return _exchange_channel( + plankton, _exchange_tracers(plankton)[$index], i, j, k, grid, fields, auxiliary_fields + ) + end +end + +# DissolvedParticulate uses `grazing` for biological removal from organic-matter pools. +@inline function grazing( + i, j, k, grid, ::Val{Tracer}, plankton::NPDPlankton, + ::NutrientsPlanktonDetritus, fields, auxiliary_fields, +) where Tracer + Tracer in _consumed_detritus(plankton) || return zero(eltype(grid)) + return -_exchange_tendency(plankton, Val(Tracer), i, j, k, grid, fields, auxiliary_fields) +end diff --git a/src/Introspection.jl b/src/Introspection.jl index 463f8801..877172fb 100644 --- a/src/Introspection.jl +++ b/src/Introspection.jl @@ -9,6 +9,7 @@ export tracer_names export auxiliary_field_names export parameter_names export parameter_domains +export model_settings export pfts export plankton_tracers export plankton_diameters @@ -21,6 +22,11 @@ export describe import Oceananigans.Biogeochemistry: required_biogeochemical_auxiliary_fields, required_biogeochemical_tracers +using ..Integrations: NPDPlankton + +@inline _introspection_target(x) = x +@inline _introspection_target(x::NPDPlankton) = x.runtime + @inline function preview_list(xs; n::Int=12) m = length(xs) @@ -40,7 +46,9 @@ underlying tracer-name tuple as a `Vector{Symbol}`. The ordering matches Oceananigans / OceanBioME state-vector conventions. """ -@inline tracer_names(bgc)::Vector{Symbol} = collect(required_biogeochemical_tracers(bgc)) +@inline tracer_names(bgc)::Vector{Symbol} = collect( + required_biogeochemical_tracers(_introspection_target(bgc)) +) """ auxiliary_field_names(bgc) -> Vector{Symbol} @@ -50,7 +58,7 @@ Auxiliary fields are non-tracer state fields (for example, light or temperature) that appear in tracer tendencies. """ @inline auxiliary_field_names(bgc)::Vector{Symbol} = collect( - required_biogeochemical_auxiliary_fields(bgc) + required_biogeochemical_auxiliary_fields(_introspection_target(bgc)) ) """ @@ -62,13 +70,24 @@ This list describes the resolved parameter fields available on the constructed biogeochemistry instance. """ function parameter_names(bgc)::Vector{Symbol} - params = getproperty(bgc, :parameters) + params = getproperty(_introspection_target(bgc), :parameters) return collect(propertynames(params)) end function _model_metadata(bgc) - hasproperty(bgc, :metadata) || return nothing - return getproperty(bgc, :metadata) + target = _introspection_target(bgc) + hasproperty(target, :metadata) || return nothing + return getproperty(target, :metadata) +end + +""" model_settings(bgc) -> NamedTuple + +Return resolved model-family properties that are not inputs to an individual process equation. +""" +function model_settings(bgc) + metadata = _model_metadata(bgc) + (metadata === nothing || !hasproperty(metadata, :model_settings)) && return NamedTuple() + return metadata.model_settings end """ pfts(bgc) -> NamedTuple @@ -152,18 +171,19 @@ function _interaction_parameter_names(bgc) end function _interaction_parameter_metadata(bgc, kind::Symbol) - available = _interaction_parameter_names(bgc) + target = _introspection_target(bgc) + available = _interaction_parameter_names(target) kind in available || begin available_text = isempty(available) ? "none" : join(string.(available), ", ") throw(ArgumentError( "Unknown interaction matrix parameter: $kind. Available parameters are: $available_text." )) end - metadata = getproperty(_model_metadata(bgc).parameter_axes, kind) - hasproperty(bgc.parameters, kind) || throw( + metadata = getproperty(_model_metadata(target).parameter_axes, kind) + hasproperty(target.parameters, kind) || throw( ArgumentError("Interaction parameter :$kind is missing from runtime parameters."), ) - matrix = getproperty(bgc.parameters, kind) + matrix = getproperty(target.parameters, kind) applicable(size, matrix) && length(size(matrix)) == 2 || throw( ArgumentError("Interaction parameter :$kind is not stored as a matrix."), ) @@ -219,11 +239,12 @@ The returned `NamedTuple` contains: - `has_sinking_velocities::Bool` """ function model_summary(bgc) + target = _introspection_target(bgc) return ( - tracers=tracer_names(bgc), - auxiliary_fields=auxiliary_field_names(bgc), - parameters=parameter_names(bgc), - has_sinking_velocities=Base.hasproperty(bgc, :sinking_velocities) && getproperty(bgc, :sinking_velocities) !== nothing, + tracers=tracer_names(target), + auxiliary_fields=auxiliary_field_names(target), + parameters=parameter_names(target), + has_sinking_velocities=Base.hasproperty(target, :sinking_velocities) && getproperty(target, :sinking_velocities) !== nothing, ) end diff --git a/src/Library/Allometry/Allometry.jl b/src/Library/Allometry/Allometry.jl index 775d94fb..b1a49635 100644 --- a/src/Library/Allometry/Allometry.jl +++ b/src/Library/Allometry/Allometry.jl @@ -2,7 +2,7 @@ module Allometry export AbstractParamDef, ConstantParam, AllometricParam -export PowerLaw +export PowerLaw, SplitPowerLaw export PalatabilityPreyParameters, PalatabilityPredatorParameters export allometric_scaling_power export allometric_palatability_unimodal, allometric_palatability_unimodal_protection diff --git a/src/Library/Allometry/interactions.jl b/src/Library/Allometry/interactions.jl index 82bdc6df..b39a84a0 100644 --- a/src/Library/Allometry/interactions.jl +++ b/src/Library/Allometry/interactions.jl @@ -167,47 +167,3 @@ function palatability_matrix_allometric_axes( return M end - -""" - consumer_assimilation_matrix_axes(T; assimilation_efficiency, - consumer_indices, prey_indices) - -Build a consumer-by-prey assimilation-efficiency matrix. - -!!! formulation - ```math - B_{ij} = \\beta_i - ``` - - where ``\\beta_i`` is the assimilation efficiency of consumer `i`. Only rows - from `consumer_indices` and columns from `prey_indices` are materialized, so - the returned matrix has size `length(consumer_indices) × length(prey_indices)`. - -# Arguments -- `T`: scalar type used for the output matrix. -- `assimilation_efficiency`: full vector of consumer assimilation efficiencies. -- `consumer_indices`: source indices for matrix rows. -- `prey_indices`: source indices for matrix columns. -""" -function consumer_assimilation_matrix_axes( - ::Type{T}; - assimilation_efficiency::AbstractVector{T}, - consumer_indices, - prey_indices, -) where {T<:Real} - n = length(assimilation_efficiency) - _validate_indices(consumer_indices, n, :consumer_indices) - _validate_indices(prey_indices, n, :prey_indices) - nr = length(consumer_indices) - nc = length(prey_indices) - M = zeros(T, nr, nc) - - @inbounds for (ii, pred) in pairs(consumer_indices) - β = assimilation_efficiency[pred] - for jj in 1:nc - M[ii, jj] = β - end - end - - return M -end diff --git a/src/Library/Allometry/parameter_defs.jl b/src/Library/Allometry/parameter_defs.jl index 83b7472b..4197aa24 100644 --- a/src/Library/Allometry/parameter_defs.jl +++ b/src/Library/Allometry/parameter_defs.jl @@ -78,6 +78,28 @@ Callable allometric power-law model using spherical cell volume. """ struct PowerLaw end +""" + SplitPowerLaw() + +Callable continuous two-regime power-law model using spherical cell volume. + +!!! formulation + ```math + p(d) = \\begin{cases} + a V(d)^{b_s}, & d \\le d_* \\\\ + a V_*^{b_s} \\left(\\frac{V(d)}{V_*}\\right)^{b_l}, & d > d_* + \\end{cases}, + \\qquad + V_* = V(d_*) + ``` + + The expected coefficient names are `prefactor` for ``a``, `breakpoint` for + the equivalent spherical diameter ``d_*``, `small_exponent` for ``b_s``, + and `large_exponent` for ``b_l``. The large-size branch is normalized at + the breakpoint so the relationship is continuous. +""" +struct SplitPowerLaw end + function allometric_relationship_identifier(model) throw( ArgumentError( @@ -87,12 +109,14 @@ function allometric_relationship_identifier(model) end allometric_relationship_identifier(::PowerLaw) = :power_law +allometric_relationship_identifier(::SplitPowerLaw) = :split_power_law function allometric_relationship_from_identifier(::Val{id}) where {id} throw(ArgumentError("Unsupported allometric relationship identifier $(repr(id)).")) end allometric_relationship_from_identifier(::Val{:power_law}) = PowerLaw() +allometric_relationship_from_identifier(::Val{:split_power_law}) = SplitPowerLaw() """ PowerLaw()(coeffs, diameter) @@ -124,6 +148,28 @@ Evaluate a `PowerLaw` allometric model. return allometric_scaling_power(a, b, diameter) end +"""Evaluate a continuous `SplitPowerLaw` at an equivalent spherical `diameter`.""" +@inline function (m::SplitPowerLaw)(coeffs::NamedTuple, diameter) + for name in (:prefactor, :breakpoint, :small_exponent, :large_exponent) + hasproperty(coeffs, name) || + throw(ArgumentError("SplitPowerLaw requires coefficient `$(name)`")) + end + + a = getproperty(coeffs, :prefactor) + breakpoint = getproperty(coeffs, :breakpoint) + small_exponent = getproperty(coeffs, :small_exponent) + large_exponent = getproperty(coeffs, :large_exponent) + breakpoint > zero(breakpoint) || + throw(ArgumentError("SplitPowerLaw `breakpoint` must be positive")) + + diameter <= breakpoint && + return allometric_scaling_power(a, small_exponent, diameter) + + value_at_breakpoint = allometric_scaling_power(a, small_exponent, breakpoint) + diameter_ratio = diameter / breakpoint + return value_at_breakpoint * diameter_ratio^(3 * large_exponent) +end + """ resolve_param(T, value, diameter) diff --git a/src/Library/nutrients.jl b/src/Library/nutrients.jl index 6a74b005..b4fe3d09 100644 --- a/src/Library/nutrients.jl +++ b/src/Library/nutrients.jl @@ -2,7 +2,7 @@ module Nutrients -export monod_limitation, liebig_minimum, frank_tnorm +export monod_limitation, inhibited_monod_limitation, liebig_minimum, frank_tnorm export normalized_droop_limitation, quota_uptake_regulation """ @@ -16,6 +16,15 @@ The indeterminate `R == K == 0` case returns zero. return R / (K + R) end +""" + inhibited_monod_limitation(resource, inhibitor, half_saturation, inhibition) + +Return a Monod resource response multiplied by exponential inhibition, +``R / (K + R) * exp(-psi I)``. +""" +@inline inhibited_monod_limitation(resource, inhibitor, half_saturation, inhibition) = + monod_limitation(resource, half_saturation) * exp(-inhibition * inhibitor) + """ liebig_minimum(a, b, rest...) liebig_minimum(values::NTuple) diff --git a/src/Library/photosynthesis.jl b/src/Library/photosynthesis.jl index 0f10e416..3b4539d8 100644 --- a/src/Library/photosynthesis.jl +++ b/src/Library/photosynthesis.jl @@ -1,7 +1,7 @@ """Light-response kernels used by phytoplankton growth formulations.""" module Photosynthesis -export smith_light_limitation, geider_light_response +export smith_light_limitation, geider_light_response, exponential_light_limitation """ smith_light_limitation(PAR, alpha, maximum_rate) @@ -23,6 +23,24 @@ slope, and `maximum_rate` is the enclosing growth-process rate scale. return light_rate / sqrt(maximum_rate * maximum_rate + light_rate * light_rate) end +""" + exponential_light_limitation(PAR, light_scale) + +Evaluate a saturating exponential light-response factor. + +```math +L_I(I) = 1 - \\exp\\left(-\\frac{I}{K_I}\\right) +``` + +`PAR` is photosynthetically active radiation and `light_scale` is the positive +irradiance scale ``K_I``. A documented ocean-biogeochemical use of this response is +Lévy, Klein & Tréguier (2001), Eq. (A7), *Journal of Marine Research* 59, 535-565, +doi:10.1357/002224001762842181. +""" +@inline exponential_light_limitation(PAR, light_scale) = + one(PAR + light_scale) - exp(-PAR / light_scale) + + """ geider_light_response(PAR, alpha, maximum_rate, chlorophyll_to_carbon_ratio) diff --git a/src/Library/predation.jl b/src/Library/predation.jl index 9f7cfc86..23f39e69 100644 --- a/src/Library/predation.jl +++ b/src/Library/predation.jl @@ -1,7 +1,7 @@ """Predation and grazing kernels.""" module Predation -export holling_type_ii, proportional_predation_loss, switching_predation_loss +export holling_type_ii, linear_predation_loss, proportional_predation_loss, switching_predation_loss """ holling_type_ii(P, K) @@ -14,6 +14,16 @@ The indeterminate `P == K == 0` case returns zero. return P / (K + P) end +""" + linear_predation_loss(inventory, consumer, grazing_rate, palatability) + +Return the prey-biomass consumption rate ``g p R Z``, where `R` is prey biomass, `Z` is +consumer biomass, `g` is the grazing coefficient, and `p` is prey palatability. Consumption +therefore increases linearly with both prey and consumer biomass. +""" +@inline linear_predation_loss(inventory, consumer, grazing_rate, palatability) = + grazing_rate * palatability * inventory * consumer + """ proportional_predation_loss( inventory, consumer, maximum_grazing_rate, half_saturation, diff --git a/src/ModelFamilies/interface.jl b/src/ModelFamilies/interface.jl index 0cc3483e..a4309576 100644 --- a/src/ModelFamilies/interface.jl +++ b/src/ModelFamilies/interface.jl @@ -1,6 +1,39 @@ export default_components export default_processes export definition_version +export plankton_roles +export ModelSetting +export setting_definitions + +"""A configurable scientific property of a model family. + +`ModelSetting` stores a default value and validation domain for model-level properties that +affect the realized model but are not parameters of an individual biological process, such as +carbon-to-nitrogen or chlorophyll-to-nitrogen ratios. +""" +struct ModelSetting{Default} + default::Default + domain::Symbol + + function ModelSetting(default::Default, domain::Symbol) where {Default} + domain in (:any, :finite, :nonnegative, :positive, :unit_interval) || throw( + ArgumentError( + "ModelSetting domain must be :any, :finite, :nonnegative, :positive, or :unit_interval" + ), + ) + return new{Default}(default, domain) + end +end + +ModelSetting(default; domain::Symbol=:finite) = ModelSetting(default, domain) + +"""Return configurable scientific properties defined by a model family.""" +setting_definitions(::AbstractModelFamily) = (;) + +"""Map user-facing plankton roles to logical components, e.g. `phytoplankton => :P`.""" +plankton_roles(::AbstractModelFamily) = throw( + ArgumentError("No method `plankton_roles(family)` is defined for this model family.") +) """Canonical logical components for a named model family. diff --git a/src/Models/FrankenLOBSTER/FrankenLOBSTER.jl b/src/Models/FrankenLOBSTER/FrankenLOBSTER.jl new file mode 100644 index 00000000..f953250f --- /dev/null +++ b/src/Models/FrankenLOBSTER/FrankenLOBSTER.jl @@ -0,0 +1,10 @@ +"""Canonical FrankenLOBSTER model family and OceanBioME plankton integration boundary.""" +module FrankenLOBSTER + +include("definition.jl") +include("parameters.jl") +include("construction.jl") + +export construct, construct_plus_recipe + +end # module diff --git a/src/Models/FrankenLOBSTER/construction.jl b/src/Models/FrankenLOBSTER/construction.jl new file mode 100644 index 00000000..0aee5658 --- /dev/null +++ b/src/Models/FrankenLOBSTER/construction.jl @@ -0,0 +1,55 @@ +using ...Construction +import ...Integrations + +function Integrations.npd_configuration(::FrankenLOBSTERFamily, _runtime) + return (; + nutrient_tracers=(:NO₃, :NH₄), consumed_detritus=(:DOM,), dependencies=(:T,), + ) +end + +function _construction_inputs(; + size_structure=DEFAULT_SIZE_STRUCTURE, + parameters::NamedTuple=(;), + settings::NamedTuple=(;), + grid=nothing, + arch=nothing, + scalar_type=nothing, + sinking_tracers=nothing, + open_bottom::Bool=true, +) + family = FrankenLOBSTERFamily() + realization = (; + plankton_pfts=Construction.plankton_realization(family, size_structure), + parameter_overrides=parameters, + setting_overrides=settings, + sinking_tracers, + open_bottom, + ) + return (; family, realization, execution=(; grid, arch, scalar_type)) +end + +"""Construct the Agate plankton component for composition with OceanBioME `LOBSTER`.""" +function construct(; kwargs...) + inputs = _construction_inputs(; kwargs...) + return Integrations.construct_npd_plankton( + inputs.family; inputs.realization..., inputs.execution... + ) +end + +"""Construct FrankenLOBSTER plankton and capture its versioned Agate recipe.""" +function construct_plus_recipe(; kwargs...) + inputs = _construction_inputs(; kwargs...) + return Integrations.construct_npd_plankton_plus_recipe( + inputs.family; inputs.realization..., inputs.execution... + ) +end + +"""Replay a FrankenLOBSTER recipe into the Agate plankton component.""" +function construct( + recipe::Construction.ModelRecipe; grid=nothing, arch=nothing, scalar_type=nothing +) + recipe.family == :FrankenLOBSTER || throw(ArgumentError( + "FrankenLOBSTER.construct requires a FrankenLOBSTER recipe; got $(recipe.family)" + )) + return Integrations.construct_npd_plankton(recipe; grid, arch, scalar_type) +end diff --git a/src/Models/FrankenLOBSTER/definition.jl b/src/Models/FrankenLOBSTER/definition.jl new file mode 100644 index 00000000..717722b2 --- /dev/null +++ b/src/Models/FrankenLOBSTER/definition.jl @@ -0,0 +1,102 @@ +using ...ModelFamilies: AbstractModelFamily +using ...Components: Plankton, Pool +using ...Processes: + Growth, Light, Consumption, Mortality, Products, ExponentialSaturation, + NutrientResponse, Monod, InhibitedMonod, Temperature, Q10, PreferentialGrazing, + HeterotrophicConsumption, LinearMortality, QuadraticMortality + +import ...ModelFamilies: default_components, default_processes, definition_version, plankton_roles +import ...Construction: family_id, registered_family + +struct FrankenLOBSTERFamily <: AbstractModelFamily end +family_id(::FrankenLOBSTERFamily) = :FrankenLOBSTER +registered_family(::Val{:FrankenLOBSTER}) = FrankenLOBSTERFamily() +definition_version(::FrankenLOBSTERFamily)::VersionNumber = v"0.14.0" +plankton_roles(::FrankenLOBSTERFamily) = ( + phytoplankton=:P, zooplankton=:Z, bacterioplankton=:H, +) + +const DEFAULT_SIZE_STRUCTURE = ( + phytoplankton=(P=(n=2, min_esd=0.6, max_esd=1.2, spacing=:linear),), + zooplankton=(Z=(n=2, min_esd=6.0, max_esd=12.0, spacing=:linear),), + bacterioplankton=(H=(n=1, min_esd=0.6, max_esd=0.6, spacing=:linear),), +) + +_nitrogen_plankton(size_structure) = Plankton(; + states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure +) + +# NO3, NH4, and DOM are external material state; T is read as a physical driver. +# Waste pools are NPD exchange accumulators. +const FRANKENLOBSTER_COMPONENTS = ( + NO₃=Pool(:nitrogen), NH₄=Pool(:nitrogen), DOM=Pool(:nitrogen), + solid_waste=Pool(:nitrogen), inorganic_waste=Pool(:nitrogen), dissolved_waste=Pool(:nitrogen), + P=_nitrogen_plankton(DEFAULT_SIZE_STRUCTURE.phytoplankton.P), + Z=_nitrogen_plankton(DEFAULT_SIZE_STRUCTURE.zooplankton.Z), + H=_nitrogen_plankton(DEFAULT_SIZE_STRUCTURE.bacterioplankton.H), +) +default_components(::FrankenLOBSTERFamily) = FRANKENLOBSTER_COMPONENTS + +const _P_GROWTH_FACTORS = ( + light=Light(ExponentialSaturation(); driver=:PAR, bindings=(light_scale=:light_half_saturation,)), + temperature=Temperature( + Q10(:plankton); driver=:T, + bindings=(q10=:temperature_q10, reference_temperature=:reference_temperature), + ), +) +const _EXUDATE_PRODUCTS = Products( + (dissolved=:dissolved_waste, inorganic=:inorganic_waste); + fractions=(inorganic=:ammonium_fraction_of_exudate,), +) + +function _p_growth(resource, response) + return Growth(; + plankton=:P, reference_resource=resource, + bindings=(maximum_rate=:maximum_growth_rate, product_fraction=:phytoplankton_exudation_fraction), + factors=merge(_P_GROWTH_FACTORS, (nutrients=response,)), products=_EXUDATE_PRODUCTS, + ) +end + +_n_mortality(plankton, rate) = Mortality( + QuadraticMortality(); plankton, bindings=(rate=rate,), products=Products(:solid_waste) +) + +const FRANKENLOBSTER_PROCESSES = ( + nitrate_growth_P=_p_growth(:NO₃, NutrientResponse( + InhibitedMonod(); resource=:NO₃, inhibitor=:NH₄, + bindings=(half_saturation=:nitrate_half_saturation, inhibition=:nitrate_ammonia_inhibition), + )), + ammonia_growth_P=_p_growth(:NH₄, NutrientResponse( + Monod(); resource=:NH₄, bindings=(half_saturation=:ammonia_half_saturation,) + )), + consumption_H_on_DOM=Consumption( + HeterotrophicConsumption(); consumers=:H, resources=:DOM, + bindings=( + maximum_rate=:bacterial_maximum_uptake_rate, + half_saturation=:bacterial_dom_half_saturation, + substrate_preference=:bacterial_substrate_preference, + assimilation=:bacterial_assimilation, + ), + unassimilated_products=:inorganic_waste, + ), + grazing_Z_on_living=Consumption( + PreferentialGrazing(); consumers=:Z, resources=(:P, :H), + bindings=( + maximum_rate=:maximum_predation_rate, half_saturation=:grazing_half_saturation, + palatability=:palatability_matrix, assimilation=:assimilation_matrix, + ), + unassimilated_products=:solid_waste, + ), + excretion_Z=Mortality( + LinearMortality(); plankton=:Z, bindings=(rate=:zooplankton_excretion_rate,), + products=Products( + (dissolved=:dissolved_waste, inorganic=:inorganic_waste); + fractions=(inorganic=:ammonium_fraction_of_zooplankton_excretion,), + ), + ), + mortality_P=_n_mortality(:P, :phytoplankton_mortality_rate), + mortality_Z=_n_mortality(:Z, :zooplankton_mortality_rate), + mortality_H=_n_mortality(:H, :bacterioplankton_mortality_rate), +) + +default_processes(::FrankenLOBSTERFamily) = FRANKENLOBSTER_PROCESSES diff --git a/src/Models/FrankenLOBSTER/parameters.jl b/src/Models/FrankenLOBSTER/parameters.jl new file mode 100644 index 00000000..209a65bb --- /dev/null +++ b/src/Models/FrankenLOBSTER/parameters.jl @@ -0,0 +1,62 @@ +import ...Parameters: + parameter_definitions, Parameter, ConstructionParameter, DerivedDefault, + DiameterIndexedVectorDefault, ConsumerResourceFromConsumer +import ...ModelFamilies: ModelSetting, setting_definitions + +using ...Library.Allometry: AllometricParam, PowerLaw +using ...Parameters: AllometricPalatability + +setting_definitions(::FrankenLOBSTERFamily) = ( + chlorophyll_ratio=ModelSetting(1.31; domain=:nonnegative), + carbon_ratio=ModelSetting(6.56; domain=:positive), +) + +"""LOBSTER3-like defaults expressed through Agate size-trait machinery.""" +function parameter_definitions(::FrankenLOBSTERFamily) + day = 86400 + law(prefactor, exponent) = AllometricParam(PowerLaw(); prefactor, exponent) + diameter_default(value) = DiameterIndexedVectorDefault(value; default=0) + + # FrankenLOBSTER keeps size-dependent traits while using LOBSTER light and N responses. + maximum_growth = law(1.2066 / day, 0.28) + nitrate_affinity = law(0.028154, 0.65) + ammonia_affinity = law(0.5 * 0.028154, 0.65) + bacterial_uptake = law(1.836 / day, 0.28) + bacterial_affinity = law(0.04284, 0.65) + + return ( + maximum_growth_rate=Parameter(diameter_default(maximum_growth)), + nitrate_half_saturation=Parameter(diameter_default(nitrate_affinity)), + ammonia_half_saturation=Parameter(diameter_default(ammonia_affinity)), + nitrate_ammonia_inhibition=Parameter(3.0), + light_half_saturation=Parameter(33.0), + temperature_q10=Parameter(1.88), + reference_temperature=Parameter(20.0), + phytoplankton_exudation_fraction=Parameter(0.05), + ammonium_fraction_of_exudate=Parameter(0.75), + phytoplankton_mortality_rate=Parameter(5.8e-7), + zooplankton_excretion_rate=Parameter(5.8e-7), + ammonium_fraction_of_zooplankton_excretion=Parameter(0.5), + zooplankton_mortality_rate=Parameter(2.31e-6), + bacterial_maximum_uptake_rate=Parameter(diameter_default(bacterial_uptake)), + bacterial_dom_half_saturation=Parameter(DerivedDefault( + ConsumerResourceFromConsumer(); deps=(:bacterial_dom_affinity_trait,) + )), + bacterial_substrate_preference=Parameter(1.0), + bacterial_assimilation=Parameter(0.1), + bacterioplankton_mortality_rate=Parameter(5.8e-7), + maximum_predation_rate=Parameter(diameter_default(law(15.9 / day, -0.16))), + grazing_half_saturation=Parameter(1.0), + palatability_matrix=Parameter(DerivedDefault( + AllometricPalatability(); deps=(:optimum_predator_prey_ratio, :specificity) + )), + assimilation_matrix=Parameter(0.7), + optimum_predator_prey_ratio=ConstructionParameter( + diameter_default(10.0); axes=:plankton + ), + specificity=ConstructionParameter(diameter_default(0.3); axes=:plankton), + bacterial_dom_affinity_trait=ConstructionParameter( + diameter_default(bacterial_affinity); axes=:plankton + ), + ) +end diff --git a/src/Models/Models.jl b/src/Models/Models.jl index afc67045..f8502033 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -6,7 +6,9 @@ module Models # ----------------------------------------------------------------------------- include("NiPiZD/NiPiZD.jl") +include("FrankenLOBSTER/FrankenLOBSTER.jl") export NiPiZD +export FrankenLOBSTER end # module diff --git a/src/Models/NiPiZD/construction.jl b/src/Models/NiPiZD/construction.jl index 89e93ce2..9b11a586 100644 --- a/src/Models/NiPiZD/construction.jl +++ b/src/Models/NiPiZD/construction.jl @@ -2,57 +2,6 @@ using OceanBioME: BoxModelGrid import ...Construction - -function _canonicalize_size_structure(size_structure) - size_structure isa NamedTuple || - throw(ArgumentError("size_structure must be a NamedTuple")) - - required_roles = (:phytoplankton, :zooplankton) - missing_roles = [role for role in required_roles if !hasproperty(size_structure, role)] - extra_roles = [role for role in keys(size_structure) if !(role in required_roles)] - isempty(missing_roles) || throw( - ArgumentError("size_structure is missing roles: $(collect(missing_roles))") - ) - isempty(extra_roles) || - throw(ArgumentError("size_structure has unknown roles: $(collect(extra_roles))")) - - phytoplankton = size_structure.phytoplankton - zooplankton = size_structure.zooplankton - phytoplankton isa NamedTuple || - throw(ArgumentError("size_structure.phytoplankton must be a NamedTuple")) - zooplankton isa NamedTuple || - throw(ArgumentError("size_structure.zooplankton must be a NamedTuple")) - isempty(phytoplankton) && - throw(ArgumentError("size_structure.phytoplankton must define at least one PFT")) - isempty(zooplankton) && - throw(ArgumentError("size_structure.zooplankton must define at least one PFT")) - - producer_pfts = keys(phytoplankton) - consumer_pfts = keys(zooplankton) - duplicate_pfts = [pft for pft in producer_pfts if pft in consumer_pfts] - isempty(duplicate_pfts) || throw( - ArgumentError( - "plankton PFT names must be unique across roles; " * - "duplicated PFTs: $(collect(duplicate_pfts))", - ), - ) - - return (; phytoplankton, zooplankton) -end - -function _plankton_realization(size_structure) - structure = _canonicalize_size_structure(size_structure) - phytoplankton = NamedTuple{keys(structure.phytoplankton)}(Tuple( - Construction.normalize_pft_size_structure(value) - for value in values(structure.phytoplankton) - )) - zooplankton = NamedTuple{keys(structure.zooplankton)}(Tuple( - Construction.normalize_pft_size_structure(value) - for value in values(structure.zooplankton) - )) - return (P=phytoplankton, Z=zooplankton) -end - function _construction_inputs(; size_structure=DEFAULT_SIZE_STRUCTURE, parameters::NamedTuple=(;), @@ -65,7 +14,7 @@ function _construction_inputs(; open_bottom::Bool=true, ) family = NiPiZDFamily() - plankton_realization = _plankton_realization(size_structure) + plankton_realization = Construction.plankton_realization(family, size_structure) parameter_overrides = parameters for (name, value) in ((:palatability_matrix, palatability_matrix), @@ -187,7 +136,7 @@ when the recipe is realized. """ function construct_plus_recipe(; kwargs...) inputs = _construction_inputs(; kwargs...) - recipe = Construction.capture_model_recipe(inputs.family; inputs.realization...) - bgc = Construction.construct(inputs.family; inputs.realization..., inputs.execution...) - return bgc, recipe + return Construction.construct_plus_recipe( + inputs.family; inputs.realization..., inputs.execution... + ) end diff --git a/src/Models/NiPiZD/definition.jl b/src/Models/NiPiZD/definition.jl index f57788b0..cdb194d9 100644 --- a/src/Models/NiPiZD/definition.jl +++ b/src/Models/NiPiZD/definition.jl @@ -5,7 +5,7 @@ using ...Processes: Smith, Monod, PreferentialGrazing, LinearMortality, QuadraticMortality, LinearRemineralization -import ...ModelFamilies: default_components, default_processes, definition_version +import ...ModelFamilies: default_components, default_processes, definition_version, plankton_roles import ...Construction: family_id, registered_family """Registered family for the size-structured NiPiZD model.""" @@ -14,6 +14,7 @@ struct NiPiZDFamily <: AbstractModelFamily end family_id(::NiPiZDFamily) = :NiPiZD registered_family(::Val{:NiPiZD}) = NiPiZDFamily() definition_version(::NiPiZDFamily)::VersionNumber = v"0.2.0" +plankton_roles(::NiPiZDFamily) = (phytoplankton=:P, zooplankton=:Z) const DEFAULT_SIZE_STRUCTURE = ( phytoplankton=(P=(n=2, min_esd=2, max_esd=10, spacing=:log),), diff --git a/src/Models/NiPiZD/parameters.jl b/src/Models/NiPiZD/parameters.jl index 3dfa16c9..07dfeb2c 100644 --- a/src/Models/NiPiZD/parameters.jl +++ b/src/Models/NiPiZD/parameters.jl @@ -13,44 +13,23 @@ import ...Parameters: using ...Library.Allometry: AllometricParam, PowerLaw -using ...Parameters: AllometricPalatability, ConsumerAssimilation +using ...Parameters: AllometricPalatability, ConsumerResourceFromConsumer function parameter_definitions(::NiPiZDFamily) detritus_remin = 0.1213 / 86400 + law(prefactor, exponent) = AllometricParam(PowerLaw(); prefactor, exponent) + diameter_default(value; default=0) = DiameterIndexedVectorDefault(value; default) return ( detritus_remineralization=Parameter(detritus_remin), mortality_export_fraction=Parameter(0.2), - linear_mortality=Parameter( - DiameterIndexedVectorDefault(8e-7; default=0) - ), - quadratic_mortality=Parameter( - DiameterIndexedVectorDefault(1e-6; default=0) - ), - maximum_growth_rate=Parameter( - DiameterIndexedVectorDefault( - AllometricParam(PowerLaw(); prefactor=2 / 86400, exponent=-0.15); - default=0, - ) - ), - nutrient_half_saturation=Parameter( - DiameterIndexedVectorDefault( - AllometricParam(PowerLaw(); prefactor=0.17, exponent=0.27); - default=0, - ) - ), - alpha=Parameter( - DiameterIndexedVectorDefault(0.1953 / 86400; default=0) - ), - maximum_predation_rate=Parameter( - DiameterIndexedVectorDefault( - AllometricParam(PowerLaw(); prefactor=30.84 / 86400, exponent=-0.16); - default=0, - ) - ), - holling_half_saturation=Parameter( - DiameterIndexedVectorDefault(5.0; default=0) - ), + linear_mortality=Parameter(diameter_default(8e-7)), + quadratic_mortality=Parameter(diameter_default(1e-6)), + maximum_growth_rate=Parameter(diameter_default(law(2 / 86400, -0.15))), + nutrient_half_saturation=Parameter(diameter_default(law(0.17, 0.27))), + alpha=Parameter(diameter_default(0.1953 / 86400)), + maximum_predation_rate=Parameter(diameter_default(law(30.84 / 86400, -0.16))), + holling_half_saturation=Parameter(diameter_default(5.0)), palatability_matrix=Parameter( DerivedDefault( AllometricPalatability(); @@ -63,24 +42,20 @@ function parameter_definitions(::NiPiZDFamily) ), assimilation_matrix=Parameter( DerivedDefault( - ConsumerAssimilation(); deps=(:assimilation_efficiency,) + ConsumerResourceFromConsumer(); deps=(:assimilation_efficiency,) ) ), optimum_predator_prey_ratio=ConstructionParameter( - DiameterIndexedVectorDefault(10.0; default=0); - axes=:plankton, + diameter_default(10.0); axes=:plankton, ), specificity=ConstructionParameter( - DiameterIndexedVectorDefault(0.3; default=0); - axes=:plankton, + diameter_default(0.3); axes=:plankton, ), protection=ConstructionParameter( - DiameterIndexedVectorDefault(0.0; default=1.0); - axes=:plankton, + diameter_default(0.0; default=1.0); axes=:plankton, ), assimilation_efficiency=ConstructionParameter( - DiameterIndexedVectorDefault(0.32; default=0); - axes=:plankton, + diameter_default(0.32); axes=:plankton, ), ) end diff --git a/src/Parameters/Parameters.jl b/src/Parameters/Parameters.jl index 4216c236..c31cf189 100644 --- a/src/Parameters/Parameters.jl +++ b/src/Parameters/Parameters.jl @@ -1,7 +1,7 @@ """Runtime parameters, construction-only parameters, and setup-time defaults.""" module Parameters -export AllometricPalatability, ConsumerAssimilation +export AllometricPalatability, ConsumerResourceFromConsumer include("parameter_types.jl") include("interaction_derivations.jl") diff --git a/src/Parameters/interaction_derivations.jl b/src/Parameters/interaction_derivations.jl index 5ef76ed6..7a1e7dc9 100644 --- a/src/Parameters/interaction_derivations.jl +++ b/src/Parameters/interaction_derivations.jl @@ -1,7 +1,6 @@ using ..Components: ModelLayout, diameter_metadata -using ..Library.Allometry: - palatability_matrix_allometric_axes, consumer_assimilation_matrix_axes +using ..Library.Allometry: palatability_matrix_allometric_axes """Return `v` when it uses the construction scalar type, otherwise throw an `ArgumentError`.""" @inline function _require_scalar_vector( @@ -16,11 +15,16 @@ using ..Library.Allometry: ) end -"""Derive consumer-by-prey palatability from allometric trait vectors.""" +"""Derive consumer-by-prey palatability from size traits, with optional prey protection.""" struct AllometricPalatability end -"""Derive consumer-by-prey assimilation from consumer-specific efficiency traits.""" -struct ConsumerAssimilation end +"""Broadcast a consumer-specific trait across a consumer-by-resource parameter matrix. + +The single declared dependency must be a vector over realized plankton SizeClasses. This is +useful when one physiological trait, such as substrate affinity, belongs to the consumer but +the runtime formulation stores a value for each consumer-resource edge. +""" +struct ConsumerResourceFromConsumer end function _plankton_entity_indices( layout::ModelLayout, labels::Tuple, parameter_name::Symbol, axis_name::Symbol @@ -52,6 +56,9 @@ end @inline function _derive_palatability(layout::ModelLayout, params, consumers, prey) _require_palatability_diameters(layout, consumers, prey) T = layout.scalar_type + protection = hasproperty(params, :protection) ? + _require_scalar_vector(T, params.protection, :protection) : + zeros(T, length(layout.size_classes)) return palatability_matrix_allometric_axes( T, layout.size_class_diameters; @@ -59,19 +66,7 @@ end T, params.optimum_predator_prey_ratio, :optimum_predator_prey_ratio ), specificity=_require_scalar_vector(T, params.specificity, :specificity), - protection=_require_scalar_vector(T, params.protection, :protection), - consumer_indices=consumers, - prey_indices=prey, - ) -end - -@inline function _derive_assimilation(layout::ModelLayout, params, consumers, prey) - T = layout.scalar_type - return consumer_assimilation_matrix_axes( - T; - assimilation_efficiency=_require_scalar_vector( - T, params.assimilation_efficiency, :assimilation_efficiency - ), + protection, consumer_indices=consumers, prey_indices=prey, ) @@ -94,17 +89,20 @@ end end @inline function _derive_parameter_default( - ::ConsumerAssimilation, + ::ConsumerResourceFromConsumer, ::Any, layout::ModelLayout, parameter, params::NamedTuple, ) + length(params) == 1 || throw(ArgumentError( + "ConsumerResourceFromConsumer requires exactly one consumer-trait dependency", + )) + trait_name = first(keys(params)) + trait = _require_scalar_vector(layout.scalar_type, first(values(params)), trait_name) consumer_labels, resource_labels = parameter.storage_labels - return _derive_assimilation( - layout, - params, - _plankton_entity_indices(layout, consumer_labels, parameter.name, :consumer), - _plankton_entity_indices(layout, resource_labels, parameter.name, :resource), + consumers = _plankton_entity_indices( + layout, consumer_labels, parameter.name, :consumer ) + return [trait[i] for i in consumers, _ in resource_labels] end diff --git a/src/Processes/Processes.jl b/src/Processes/Processes.jl index 67eecec9..c6b09d54 100644 --- a/src/Processes/Processes.jl +++ b/src/Processes/Processes.jl @@ -6,17 +6,17 @@ using ..Components: Plankton, Pool, PlanktonStateRef, ModelLayout, element, stat using ..ModelFamilies: AbstractModelFamily, default_components, default_processes using ..Parameters: Parameter, ConstructionParameter, DerivedDefault, parameter_definitions using ..Library.Mortality: linear_loss -using ..Library.Predation: proportional_predation_loss, switching_predation_loss -using ..Library.Photosynthesis: geider_light_response, smith_light_limitation +using ..Library.Predation: linear_predation_loss, proportional_predation_loss, switching_predation_loss +using ..Library.Photosynthesis: exponential_light_limitation, geider_light_response, smith_light_limitation using ..Library.Nutrients: - frank_tnorm, liebig_minimum, monod_limitation, normalized_droop_limitation, - quota_uptake_regulation + frank_tnorm, inhibited_monod_limitation, liebig_minimum, monod_limitation, + normalized_droop_limitation, quota_uptake_regulation using ..Library.Temperature: q10_temperature_factor using ..Library.Remineralization: linear_remineralization export AbstractProcess, AbstractFormulation, AbstractFactor, AbstractStoichiometry -export Smith, Geider, Monod, NormalizedDroop, QuotaRegulatedMonod, Liebig, FrankTNorm, Q10 -export PreferentialGrazing, HeterotrophicConsumption +export Smith, Geider, ExponentialSaturation, Monod, InhibitedMonod, NormalizedDroop, QuotaRegulatedMonod, Liebig, FrankTNorm, Q10 +export PreferentialGrazing, LinearGrazing, HeterotrophicConsumption export LinearMortality, QuadraticMortality, LinearRemineralization export Light, NutrientLimitation, Temperature export FixedStoichiometry diff --git a/src/Processes/canonical_semantics.jl b/src/Processes/canonical_semantics.jl index f5f109b3..0e339c4d 100644 --- a/src/Processes/canonical_semantics.jl +++ b/src/Processes/canonical_semantics.jl @@ -165,7 +165,42 @@ function process_facts(process::Growth, id::Symbol, components::NamedTuple) end end - return (; plankton_states) + product_targets = if isnothing(process.products) + nothing + else + targets = _canonical_product_targets( + id, process.products, components, reference_element, "growth products" + ) + _product_transfer_mode( + id, + process.products, + targets, + (reference_element,), + reference_element, + "growth products", + ) + if has_stoichiometry + product_stoichiometry = process.products.stoichiometry + isnothing(product_stoichiometry) && throw(ArgumentError( + "process :$id fixed-stoichiometry growth products must use FixedStoichiometry " * + "so routed products account for every growth element", + )) + growth_elements = (reference_element, keys(process.additional_resources)...) + product_elements = Tuple(keys(first(values(targets)))) + sort(collect(product_elements); by=String) == sort(collect(growth_elements); by=String) || + throw(ArgumentError( + "process :$id fixed-stoichiometry growth products must route elements " * + "$growth_elements; got $product_elements", + )) + product_stoichiometry.bindings == process.stoichiometry.bindings || throw(ArgumentError( + "process :$id growth-product stoichiometry must use the same ratio bindings as " * + "growth stoichiometry", + )) + end + targets + end + + return (; plankton_states, product_targets) end function process_facts( diff --git a/src/Processes/factor_vocabulary.jl b/src/Processes/factor_vocabulary.jl index 0709b7fa..9e6db5d6 100644 --- a/src/Processes/factor_vocabulary.jl +++ b/src/Processes/factor_vocabulary.jl @@ -14,9 +14,15 @@ synthesize a prognostic non-elemental state such as `:chlorophyll`. """ struct Geider <: AbstractFormulation end -"""Monod single-resource limitation formulation.""" +"""Saturating-exponential light limitation, ``1 - exp(-I / K_I)``.""" +struct ExponentialSaturation <: AbstractFormulation end + +"""Monod saturation formulation, ``x / (K + x)``.""" struct Monod <: AbstractFormulation end +"""Monod resource limitation multiplied by exponential inhibition.""" +struct InhibitedMonod <: AbstractFormulation end + """Normalized Droop cellular-quota growth-limitation formulation.""" struct NormalizedDroop <: AbstractFormulation end @@ -29,8 +35,9 @@ struct Liebig <: AbstractFormulation end """Differentiable Frank t-norm nutrient-combination formulation.""" struct FrankTNorm <: AbstractFormulation end -"""Q10 temperature-response formulation.""" -struct Q10 <: AbstractFormulation end +"""Q10 temperature-response formulation indexed over the affected process participant role.""" +struct Q10{Axis} <: AbstractFormulation end +Q10(axis::Symbol) = Q10{axis}() """Growth formulation with a base maximum rate and optional multiplicative factors.""" struct FactorizedGrowth <: AbstractFormulation end @@ -64,10 +71,18 @@ function PreferentialGrazing(; switching_exponent=1) return PreferentialGrazing(switching_exponent) end +"""Linear grazing with prey consumption proportional to prey and consumer biomass. + +The `rate` parameter is consumer-indexed and has inverse-concentration inverse-time units. +`palatability` scales individual consumer-resource links. +""" +struct LinearGrazing <: AbstractFormulation end + """Heterotrophic consumption of substitutable substrates with shared consumer capacity. `maximum_rate` is one per-consumer uptake capacity shared across all declared substrates. -`substrate_preference` controls the relative accessibility of each consumer-resource pair. +`half_saturation` and `substrate_preference` are consumer-resource properties, allowing consumers +to differ in affinity and relative accessibility for the same substrate. """ struct HeterotrophicConsumption <: AbstractFormulation end @@ -142,14 +157,15 @@ function _canonical_participants(role::Symbol, values) end """Light-dependent multiplicative Growth factor using the Growth rate scale.""" -struct Light{Formulation<:Union{Smith,Geider}} <: AbstractFactor +struct Light{Formulation<:Union{Smith,Geider,ExponentialSaturation,Monod}} <: AbstractFactor formulation::Formulation driver::Symbol bindings::NamedTuple end function Light( - formulation::Union{Smith,Geider}; driver::Symbol, bindings::NamedTuple=NamedTuple() + formulation::Union{Smith,Geider,ExponentialSaturation,Monod}; + driver::Symbol, bindings::NamedTuple=NamedTuple(), ) return Light(formulation, driver, _canonical_bindings(bindings)) end @@ -160,16 +176,26 @@ authored_parameter_bindings(factor::Light) = factor.bindings The factor reads an environmental Pool but does not define process material transfer. """ -struct NutrientResponse{Formulation<:Monod} <: AbstractFactor +struct NutrientResponse{Formulation<:Union{Monod,InhibitedMonod},Inhibitor} <: AbstractFactor formulation::Formulation resource::Symbol + inhibitor::Inhibitor bindings::NamedTuple end function NutrientResponse( formulation::Monod; resource::Symbol, bindings::NamedTuple=NamedTuple() ) - return NutrientResponse(formulation, resource, _canonical_bindings(bindings)) + return NutrientResponse(formulation, resource, nothing, _canonical_bindings(bindings)) +end + +function NutrientResponse( + formulation::InhibitedMonod; + resource::Symbol, inhibitor::Symbol, bindings::NamedTuple=NamedTuple(), +) + return NutrientResponse( + formulation, resource, inhibitor, _canonical_bindings(bindings) + ) end authored_parameter_bindings(factor::NutrientResponse) = factor.bindings @@ -195,17 +221,29 @@ end authored_parameter_bindings(factor::QuotaResponse) = factor.bindings -"""Temperature-dependent multiplicative process-rate factor.""" -struct Temperature{Formulation<:Q10} <: AbstractFactor +"""Temperature-dependent multiplicative process-rate factor. + +By default temperature is read from an external driver named `:temperature`. Pass `component` +instead to read a scalar model component such as an Oceananigans temperature tracer. +""" +struct Temperature{Formulation<:Q10,Driver,Component} <: AbstractFactor formulation::Formulation - driver::Symbol + driver::Driver + component::Component bindings::NamedTuple end function Temperature( - formulation::Q10; driver::Symbol=:temperature, bindings::NamedTuple=NamedTuple() + formulation::Q10; + driver::Union{Nothing,Symbol}=nothing, + component::Union{Nothing,Symbol}=nothing, + bindings::NamedTuple=NamedTuple(), ) - return Temperature(formulation, driver, _canonical_bindings(bindings)) + isnothing(driver) || isnothing(component) || throw( + ArgumentError("Temperature accepts either `driver` or `component`, not both"), + ) + isnothing(driver) && isnothing(component) && (driver = :temperature) + return Temperature(formulation, driver, component, _canonical_bindings(bindings)) end authored_parameter_bindings(factor::Temperature) = factor.bindings @@ -278,8 +316,12 @@ end """Return the ordered semantic inputs read by a factor before its parameter slots.""" factor_inputs(::AbstractFactor) = () factor_inputs(factor::Light) = (FactorDriver(factor.driver),) -factor_inputs(factor::Temperature) = (FactorDriver(factor.driver),) -factor_inputs(factor::NutrientResponse) = (FactorComponent(factor.resource),) +factor_inputs(factor::Temperature) = isnothing(factor.component) ? + (FactorDriver(factor.driver),) : (FactorComponent(factor.component),) +factor_inputs(factor::NutrientResponse{<:Monod}) = (FactorComponent(factor.resource),) +factor_inputs(factor::NutrientResponse{<:InhibitedMonod}) = ( + FactorComponent(factor.resource), FactorComponent(factor.inhibitor), +) factor_inputs(::QuotaResponse) = () """Return named child factors composed by a factor.""" diff --git a/src/Processes/parameter_schema.jl b/src/Processes/parameter_schema.jl index 22300e88..6987576f 100644 --- a/src/Processes/parameter_schema.jl +++ b/src/Processes/parameter_schema.jl @@ -45,14 +45,26 @@ parameter_slots(::AbstractFormulation) = () parameter_slots(::FactorizedGrowth) = ( ParameterSlot(:maximum_rate, (:plankton,); domain=:nonnegative), ) +parameter_slots(process::Growth) = isnothing(process.products) ? + parameter_slots(FactorizedGrowth()) : ( + parameter_slots(FactorizedGrowth())..., + ParameterSlot(:product_fraction, (:plankton,); domain=:unit_interval), + ) parameter_slots(::Smith) = (ParameterSlot(:alpha, (:plankton,); domain=:nonnegative),) parameter_slots(::Geider) = ( ParameterSlot(:alpha, (:plankton,); domain=:nonnegative), ParameterSlot(:chlorophyll_to_carbon_ratio, (:plankton,); domain=:nonnegative), ) +parameter_slots(::ExponentialSaturation) = ( + ParameterSlot(:light_scale, (:plankton,); domain=:positive), +) parameter_slots(::Monod) = ( ParameterSlot(:half_saturation, (:plankton,); domain=:nonnegative), ) +parameter_slots(::InhibitedMonod) = ( + ParameterSlot(:half_saturation, (:plankton,); domain=:nonnegative), + ParameterSlot(:inhibition; domain=:nonnegative), +) parameter_slots(::NormalizedDroop) = ( ParameterSlot(:minimum_quota, (:plankton,); domain=:positive), ParameterSlot(:maximum_quota, (:plankton,); domain=:positive), @@ -66,8 +78,8 @@ parameter_slots(::QuotaRegulatedMonod) = ( ) parameter_slots(::Liebig) = () parameter_slots(::FrankTNorm) = (ParameterSlot(:sharpness; domain=:positive),) -parameter_slots(::Q10) = ( - ParameterSlot(:q10; domain=:positive), +parameter_slots(::Q10{Axis}) where {Axis} = ( + ParameterSlot(:q10, (Axis,); domain=:positive), ParameterSlot(:reference_temperature), ) parameter_slots(::PreferentialGrazing) = ( @@ -76,9 +88,14 @@ parameter_slots(::PreferentialGrazing) = ( ParameterSlot(:palatability, (:consumer, :resource); domain=:nonnegative), ParameterSlot(:assimilation, (:consumer, :resource); domain=:unit_interval), ) +parameter_slots(::LinearGrazing) = ( + ParameterSlot(:rate, (:consumer,); domain=:nonnegative), + ParameterSlot(:palatability, (:consumer, :resource); domain=:nonnegative), + ParameterSlot(:assimilation, (:consumer, :resource); domain=:unit_interval), +) parameter_slots(::HeterotrophicConsumption) = ( ParameterSlot(:maximum_rate, (:consumer,); domain=:nonnegative), - ParameterSlot(:half_saturation, (:resource,); domain=:positive), + ParameterSlot(:half_saturation, (:consumer, :resource); domain=:positive), ParameterSlot(:substrate_preference, (:consumer, :resource); domain=:nonnegative), ParameterSlot(:assimilation, (:consumer, :resource); domain=:unit_interval), ) @@ -113,6 +130,8 @@ struct ParameterBinding{Axes,AxisComponents} domain::Symbol end -_parameter_slot_source(node::Union{AbstractFormulation,AbstractStoichiometry,Products}) = node +_parameter_slot_source( + node::Union{AbstractFormulation,AbstractStoichiometry,Products,Growth} +) = node _parameter_slot_source(node) = formulation(node) diff --git a/src/Processes/parameter_validation.jl b/src/Processes/parameter_validation.jl index 7c5becc7..35009815 100644 --- a/src/Processes/parameter_validation.jl +++ b/src/Processes/parameter_validation.jl @@ -1,4 +1,5 @@ @inline function parameter_domain_valid(value, domain::Symbol) + domain === :any && return true value isa Real && !(value isa Bool) && isfinite(value) || return false domain === :finite && return true domain === :nonnegative && return value >= zero(value) diff --git a/src/Processes/process_declarations.jl b/src/Processes/process_declarations.jl index c1b83238..c909c362 100644 --- a/src/Processes/process_declarations.jl +++ b/src/Processes/process_declarations.jl @@ -117,15 +117,21 @@ end `bindings.maximum_rate` names the model parameter that sets the growth-rate scale. `reference_resource` supplies the Element represented by the plankton `reference_state`. `additional_resources` maps additional Elements to external Pools consumed according to -`FixedStoichiometry`. Factors modify growth rate only; independently prognostic elemental -states are supplied through [`NutrientUptake`](@ref). +`FixedStoichiometry`. Factors modify gross growth rate only. Optional `products` route the +`product_fraction` of gross growth before biomass retention while resource uptake remains gross; +the retained biomass fraction is the exact complement. For fixed-stoichiometry growth, routed +products must account for every growth element using the same stoichiometric ratio bindings. +Independently prognostic elemental states are supplied through [`NutrientUptake`](@ref). """ -struct Growth{Factors<:NamedTuple,AdditionalResources<:NamedTuple,Stoichiometry} <: AbstractProcess +struct Growth{ + Factors<:NamedTuple,AdditionalResources<:NamedTuple,Stoichiometry,ProductRouting +} <: AbstractProcess plankton::Tuple factors::Factors reference_resource::Symbol additional_resources::AdditionalResources stoichiometry::Stoichiometry + products::ProductRouting bindings::NamedTuple end @@ -135,6 +141,7 @@ function Growth(; factors::NamedTuple=NamedTuple(), additional_resources::NamedTuple=NamedTuple(), stoichiometry=nothing, + products=nothing, bindings::NamedTuple=NamedTuple(), ) all(resource -> resource isa Symbol, values(additional_resources)) || throw( @@ -149,6 +156,7 @@ function Growth(; reference_resource, _canonical_namedtuple(additional_resources), stoichiometry, + _canonical_products(products), _canonical_bindings(bindings), ) end @@ -190,13 +198,14 @@ authored_parameter_bindings(process::NutrientUptake) = process.bindings """Consumer-resource process with optional factors and unassimilated products. For `PreferentialGrazing`, `maximum_rate` is one consumer-level ingestion capacity shared across -all declared prey. For `HeterotrophicConsumption`, `maximum_rate` is likewise one consumer-level -uptake capacity shared across substitutable substrates. When one living-prey consumption process -routes multi-element unassimilated products from multiple resources, those resources currently +all declared prey. `LinearGrazing` instead applies a consumer-level mass-action `rate` +independently to each consumer-resource link. For `HeterotrophicConsumption`, `maximum_rate` is +likewise one consumer-level uptake capacity shared across substitutable substrates. When one living-prey +consumption process routes multi-element unassimilated products from multiple resources, those resources currently must expose the same prognostic Element set. """ struct Consumption{ - Formulation<:Union{PreferentialGrazing,HeterotrophicConsumption}, + Formulation<:Union{PreferentialGrazing,LinearGrazing,HeterotrophicConsumption}, Factors<:NamedTuple, ProductRouting, } <: AbstractProcess @@ -209,7 +218,7 @@ struct Consumption{ end function Consumption( - formulation::Union{PreferentialGrazing,HeterotrophicConsumption}; + formulation::Union{PreferentialGrazing,LinearGrazing,HeterotrophicConsumption}; consumers, resources, factors::NamedTuple=NamedTuple(), @@ -289,13 +298,14 @@ factors(::AbstractProcess) = NamedTuple() factors(process::Union{Growth,Consumption}) = process.factors process_products(::AbstractProcess) = nothing -process_products(process::Union{Consumption,Mortality}) = process.products -product_path(::Mortality) = (:products,) +process_products(process::Union{Growth,Consumption,Mortality}) = process.products +product_path(::Union{Growth,Mortality}) = (:products,) product_path(::Consumption) = (:unassimilated_products,) """Whether a consumer-resource formulation uses living consumer-prey interaction matrices.""" uses_living_interactions(::AbstractFormulation) = false uses_living_interactions(::PreferentialGrazing) = true +uses_living_interactions(::LinearGrazing) = true """Return canonical participant roles for an authored scientific process.""" function participants(process::Growth) diff --git a/src/Processes/rates.jl b/src/Processes/rates.jl index 54c91749..deb3e6c7 100644 --- a/src/Processes/rates.jl +++ b/src/Processes/rates.jl @@ -18,9 +18,16 @@ function factor_value end ::Geider, light, maximum_rate, alpha, chlorophyll_to_carbon_ratio ) = geider_light_response(light, alpha, maximum_rate, chlorophyll_to_carbon_ratio) +@inline factor_value(::ExponentialSaturation, light, light_scale) = + exponential_light_limitation(light, light_scale) + @inline factor_value(::Monod, resource, half_saturation) = monod_limitation(resource, half_saturation) +@inline factor_value( + ::InhibitedMonod, resource, inhibitor, half_saturation, inhibition +) = inhibited_monod_limitation(resource, inhibitor, half_saturation, inhibition) + @inline factor_value( ::NormalizedDroop, internal, reference, minimum_quota, maximum_quota ) = normalized_droop_limitation(internal, reference, minimum_quota, maximum_quota) @@ -54,6 +61,11 @@ function factor_value end ) end +"""Evaluate mass-action loss of one living prey state.""" +@inline process_rate( + ::LinearGrazing, inventory, consumer, rate, palatability +) = linear_predation_loss(inventory, consumer, rate, palatability) + """Evaluate one substrate uptake rate from a shared heterotrophic consumer capacity.""" @inline function process_rate( ::HeterotrophicConsumption, diff --git a/test/runtests.jl b/test/runtests.jl index 86e5204d..6ad116c0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ include("test_multistate_process_compilation.jl") include("test_food_web_compilation.jl") include("test_direct_model_definition_construction.jl") include("test_models_construct.jl") +include("test_frankenlobster.jl") include("test_forwarddiff.jl") include("test_active_parameters.jl") include("test_enzyme.jl") diff --git a/test/test_direct_model_definition_construction.jl b/test/test_direct_model_definition_construction.jl index a852e267..4fab4df6 100644 --- a/test/test_direct_model_definition_construction.jl +++ b/test/test_direct_model_definition_construction.jl @@ -3,7 +3,7 @@ using Oceananigans.Biogeochemistry: using Test using Agate.Components: Plankton, Pool -using Agate.Parameters: AllometricPalatability, ConsumerAssimilation, ConstantDefault, ConstructionParameter, DerivedDefault, Parameter +using Agate.Parameters: AllometricPalatability, ConsumerResourceFromConsumer, ConstantDefault, ConstructionParameter, DerivedDefault, Parameter using Agate.Construction: construct using Agate.Introspection: plankton_diameters using Agate.Processes: @@ -71,7 +71,7 @@ function direct_npz_definition() ) ), assimilation_matrix=Parameter( - DerivedDefault(ConsumerAssimilation(); deps=(:assimilation_efficiency,)) + DerivedDefault(ConsumerResourceFromConsumer(); deps=(:assimilation_efficiency,)) ), ) return ModelDefinition(; components, processes, parameters) @@ -114,7 +114,7 @@ end parameters = ( maximum_predation_rate=Parameter(1.0), half_saturation=Parameter(1.0), palatability_shared=Parameter(0.5), assimilation_other=Parameter(0.5), - assimilation_local=Parameter(DerivedDefault(ConsumerAssimilation(); deps=(:assimilation_efficiency,))), + assimilation_local=Parameter(DerivedDefault(ConsumerResourceFromConsumer(); deps=(:assimilation_efficiency,))), assimilation_efficiency=ConstructionParameter(0.5; axes=:plankton), ) bgc = construct( diff --git a/test/test_food_web_compilation.jl b/test/test_food_web_compilation.jl index 9e7e5ddc..6afb37bf 100644 --- a/test/test_food_web_compilation.jl +++ b/test/test_food_web_compilation.jl @@ -7,7 +7,7 @@ using Agate.Construction: construct using Agate.Parameters: Parameter, NoDefault using Agate.Processes: ModelDefinition, Growth, Light, NutrientResponse, Temperature, Consumption, Smith, Monod, - Q10, HeterotrophicConsumption, PreferentialGrazing, participants + Q10, HeterotrophicConsumption, PreferentialGrazing, LinearGrazing, participants function food_web_definition(; grazing=PreferentialGrazing()) components = ( @@ -19,8 +19,13 @@ function food_web_definition(; grazing=PreferentialGrazing()) M=Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[2.0]), Z=Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[10.0]), ) - temperature = Temperature( - Q10(); bindings=(q10=:temperature_q10, reference_temperature=:reference_temperature) + plankton_temperature = Temperature( + Q10(:plankton); + bindings=(q10=:plankton_temperature_q10, reference_temperature=:reference_temperature), + ) + consumer_temperature = Temperature( + Q10(:consumer); + bindings=(q10=:consumer_temperature_q10, reference_temperature=:reference_temperature), ) processes = ( growth_autotrophs=Growth(; @@ -28,7 +33,7 @@ function food_web_definition(; grazing=PreferentialGrazing()) reference_resource=:N, bindings=(maximum_rate=:maximum_growth_rate,), factors=( - temperature=temperature, + temperature=plankton_temperature, nutrients=NutrientResponse( Monod(); resource=:N, bindings=(half_saturation=:nutrient_half_saturation,) @@ -46,7 +51,7 @@ function food_web_definition(; grazing=PreferentialGrazing()) substrate_preference=:substrate_preference_matrix, assimilation=:bacterial_assimilation, ), - factors=(temperature=temperature,), + factors=(temperature=consumer_temperature,), unassimilated_products=:D, ), grazing_living=Consumption( @@ -67,7 +72,8 @@ function food_web_definition(; grazing=PreferentialGrazing()) maximum_growth_rate=no_default(), alpha=no_default(), nutrient_half_saturation=no_default(), - temperature_q10=no_default(), + plankton_temperature_q10=no_default(), + consumer_temperature_q10=no_default(), reference_temperature=no_default(), maximum_consumption_rate=no_default(), pom_half_saturation=no_default(), @@ -86,10 +92,11 @@ function food_web_parameter_overrides(::Type{T}=Float64) where {T<:Real} maximum_growth_rate=T[2e-5, 1.4e-5], alpha=T[2e-6, 1.6e-6], nutrient_half_saturation=T[0.2, 0.3], - temperature_q10=T(2), + plankton_temperature_q10=T[2, 2], + consumer_temperature_q10=T[2], reference_temperature=T(20), maximum_consumption_rate=T[1.5e-5], - pom_half_saturation=T[0.15], + pom_half_saturation=reshape(T[0.15], 1, 1), substrate_preference_matrix=reshape(T[1.0], 1, 1), bacterial_assimilation=reshape(T[0.65], 1, 1), maximum_predation_rate=T[6e-5, 9e-5], @@ -120,7 +127,8 @@ end (:maximum_growth_rate, [NaN, 1.4e-5], :nonnegative, "NaN"), (:maximum_growth_rate, [-1.0, 1.4e-5], :nonnegative, "-1.0"), (:reference_temperature, Inf, :finite, "Inf"), - (:temperature_q10, 0.0, :positive, "0.0"), + (:plankton_temperature_q10, [0.0, 2.0], :positive, "0.0"), + (:consumer_temperature_q10, [0.0], :positive, "0.0"), (:living_palatability_matrix, [NaN 0.8; 0.7 0.9], :nonnegative, "NaN"), (:living_assimilation_matrix, [-0.1 0.5; 0.35 0.45], :unit_interval, "-0.1"), (:living_assimilation_matrix, [1.1 0.5; 0.35 0.45], :unit_interval, "1.1"), @@ -171,7 +179,7 @@ end ) @test process_compiler_isapprox(growth30, 2 * growth20) direct_growth20 = 0.05 * 2e-5 * - Agate.Processes.factor_value(Q10(), 20.0, 2.0, 20.0) * + Agate.Processes.factor_value(Q10(:plankton), 20.0, 2.0, 20.0) * Agate.Processes.factor_value(Monod(), 5.0, 0.2) * Agate.Processes.factor_value(Smith(), 100.0, 2e-5, 2e-6) @test process_compiler_isapprox(growth20, direct_growth20) @@ -215,7 +223,7 @@ end definition = ModelDefinition(; components, processes, parameters) base_overrides = ( maximum_consumption_rate=[2.0], - pom_half_saturation=[1.0, 3.0, 7.0], + pom_half_saturation=reshape([1.0, 3.0, 7.0], 1, 3), substrate_preference_matrix=ones(1, 3), bacterial_assimilation=reshape([0.2, 0.4, 0.8], 1, 3), ) @@ -281,3 +289,31 @@ end @test prey_losses(zero_proportional, 0.0, 0.0) == (0.0, 0.0) @test prey_losses(zero_switching, 0.0, 0.0) == (0.0, 0.0) end + +@testset "Linear grazing is exact mass action" begin + components = ( + P=Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[1.0]), + B=Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[0.8]), + Z=Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[10.0]), + D=Pool(:nitrogen), + ) + processes = (grazing=Consumption( + LinearGrazing(); consumers=:Z, resources=(:P, :B), + bindings=(rate=:grazing_rate, palatability=:palatability, assimilation=:assimilation), + unassimilated_products=:D, + ),) + parameters = ( + grazing_rate=Parameter(0.25), + palatability=Parameter([0.8 0.5]), + assimilation=Parameter([0.7 0.6]), + ) + model = construct(Agate.Processes.ModelDefinition(; components, processes, parameters)) + args = food_web_args(model, (P_1=2.0, B_1=4.0, Z_1=3.0, D=0.0)) + p_loss, b_loss = 0.25 * 0.8 * 2.0 * 3.0, 0.25 * 0.5 * 4.0 * 3.0 + actual = ( + model(Val(:P_1), args...), model(Val(:B_1), args...), + model(Val(:Z_1), args...), model(Val(:D), args...), + ) + expected = (-p_loss, -b_loss, 0.7 * p_loss + 0.6 * b_loss, 0.3 * p_loss + 0.4 * b_loss) + @test all(isapprox.(actual, expected)) +end diff --git a/test/test_frankenlobster.jl b/test/test_frankenlobster.jl new file mode 100644 index 00000000..d8750447 --- /dev/null +++ b/test/test_frankenlobster.jl @@ -0,0 +1,134 @@ +using Test +using Oceananigans.Architectures: CPU +using Oceananigans.Grids: RectilinearGrid +using Oceananigans.Fields: ConstantField +using Oceananigans.Biogeochemistry: + required_biogeochemical_auxiliary_fields, required_biogeochemical_tracers + +using OceanBioME: chlorophyll, PrescribedPhotosyntheticallyActiveRadiation +using OceanBioME.Models.NutrientsPlanktonDetritusModels: + DissolvedParticulate, LOBSTER, nutrient_uptake +using OceanBioME.Models.NutrientsPlanktonDetritusModels.NutrientsModels: + Nutrients, NitrateAmmonia + +const FrankenLOBSTER = Agate.Models.FrankenLOBSTER +const _GRID = RectilinearGrid(CPU(); size=(1, 1, 1), extent=(1, 1, 1)) +_cell(x) = fill(x, 1, 1, 1) +_light(x=1.0) = PrescribedPhotosyntheticallyActiveRadiation(ConstantField(x)) + +function _fields(; NO₃=0.0, NH₄=0.0, T=20.0, DOM=0.0, sPOM=0.0, bPOM=0.0, + P_1=0.0, P_2=0.0, Z_1=0.0, Z_2=0.0, H_1=0.0) + state = (; NO₃, NH₄, T, DOM, sPOM, bPOM, P_1, P_2, Z_1, Z_2, H_1) + return NamedTuple{keys(state)}(map(_cell, values(state))) +end + +const _CONTROLLED = ( + maximum_growth_rate=(P_1=1.0, P_2=1.0), nitrate_half_saturation=(P_1=1.0, P_2=1.0), + ammonia_half_saturation=(P_1=1.0, P_2=1.0), nitrate_ammonia_inhibition=0.1, + light_half_saturation=(P_1=1.0, P_2=1.0), + temperature_q10=(P_1=2.0, P_2=2.0), reference_temperature=20.0, + phytoplankton_mortality_rate=(P_1=0.0, P_2=0.0), maximum_predation_rate=(Z_1=0.0, Z_2=0.0), + zooplankton_excretion_rate=(Z_1=1.0, Z_2=1.0), zooplankton_mortality_rate=(Z_1=0.0, Z_2=0.0), + bacterial_maximum_uptake_rate=(H_1=2.0,), bacterial_dom_half_saturation=reshape([1.0], 1, 1), + bacterial_substrate_preference=reshape([1.0], 1, 1), bacterial_assimilation=reshape([0.25], 1, 1), + bacterioplankton_mortality_rate=(H_1=0.0,), +) + +function _controlled(; parameters=(;)) + plankton = FrankenLOBSTER.construct(; grid=_GRID, parameters=merge(_CONTROLLED, parameters)) + detritus = DissolvedParticulate( + _GRID; dissolved_remineralisation_rate=0.0, + particulate_remineralisation_rate=(0.0, 0.0), sinking_speeds=(0.0, 0.0), + ) + return LOBSTER( + _GRID; plankton, nutrients=Nutrients(; nitrogen=NitrateAmmonia(; nitrification_rate=0.0)), + light_attenuation=_light(), detritus, + ).underlying_biogeochemistry +end + +@testset "FrankenLOBSTER composition and replay" begin + size_structure = ( + phytoplankton=(pico=[0.5], nano=[2.0]), zooplankton=(micro=[8.0], meso=[20.0]), + bacterioplankton=(heterotroph=[0.4, 0.8],), + ) + parameters = ( + assimilation_matrix=fill(0.65, 2, 4), maximum_growth_rate=(nano_1=1e-5,), + ) + settings = (chlorophyll_ratio=1.5,) + plankton, recipe = FrankenLOBSTER.construct_plus_recipe(; + grid=_GRID, arch=CPU(), scalar_type=Float32, size_structure, parameters, settings, + sinking_tracers=(nano_1=0.1,), open_bottom=false, + ) + decoded = Agate.Construction.decode_recipe(Agate.Construction.encode_recipe(recipe)) + replayed = FrankenLOBSTER.construct(decoded; grid=_GRID, arch=CPU(), scalar_type=Float32) + direct32 = FrankenLOBSTER.construct(; arch=CPU(), scalar_type=Float32) + bgc = LOBSTER(_GRID; plankton) + + @test required_biogeochemical_tracers(plankton) == + (:nano_1, :pico_1, :meso_1, :micro_1, :heterotroph_1, :heterotroph_2) + @test size(plankton.runtime.parameters.bacterial_dom_half_saturation) == (2, 1) + @test :T ∉ required_biogeochemical_tracers(plankton.runtime) + @test required_biogeochemical_auxiliary_fields(plankton.runtime) == (:PAR, :T) + @test required_biogeochemical_auxiliary_fields(plankton) == (:PAR,) + @test all(t -> t in required_biogeochemical_tracers(bgc), (:NO₃, :NH₄, :DOM, :sPOM, :bPOM, :T)) + @test chlorophyll(plankton, (tracers=(nano_1=_cell(2.0), pico_1=_cell(1.0)),))[1, 1, 1] ≈ 4.5 + @test recipe.setting_overrides == settings + + I = Agate.Introspection + for inspect in ( + I.model_settings, I.pfts, I.parameter_names, I.plankton_tracers, + I.plankton_diameters, I.model_summary, + ) + @test inspect(plankton) == inspect(plankton.runtime) + end + @test I.parameter_domains(plankton, :assimilation_matrix) == + I.parameter_domains(plankton.runtime, :assimilation_matrix) + @test I.interaction_matrix(plankton, :assimilation_matrix) == + I.interaction_matrix(plankton.runtime, :assimilation_matrix) + + @test (replayed.runtime.parameters, replayed.traits) == (plankton.runtime.parameters, plankton.traits) + @test eltype(plankton.runtime.parameters.maximum_growth_rate) === Float32 + @test eltype(direct32.runtime.parameters.maximum_growth_rate) === Float32 + @test_throws ArgumentError Agate.Integrations.NPDPlankton( + plankton.runtime; + owned_components=(:P, :Z, :H), + phytoplankton_components=(:P,), + exchange_tracers=(solid=:solid_watse, dissolved=nothing, inorganic=nothing), + traits=plankton.traits, + ) + @test Agate.Integrations.NPDPlankton( + plankton.runtime; + owned_components=(:P, :Z, :H), + phytoplankton_components=(:P,), + nutrient_tracers=(:DOM,), + traits=plankton.traits, + ) isa Agate.Integrations.NPDPlankton + @test_throws ArgumentError FrankenLOBSTER.construct(settings=(unknown=1.0,)) + @test_throws ArgumentError FrankenLOBSTER.construct(sinking_tracers=(P_1=0.1,)) +end + +@testset "FrankenLOBSTER LOBSTER exchange" begin + bgc = _controlled() + aux = (PAR=_cell(1.0),) + tendency(tracer, fields) = bgc(1, 1, 1, _GRID, Val(tracer), (; time=0.0), fields, aux) + uptake(tracer, fields) = nutrient_uptake(1, 1, 1, _GRID, Val(tracer), bgc.plankton, bgc, fields, aux) + + @test uptake(:DOM, _fields()) == 0 + + nitrate = _fields(; NO₃=1.0, P_1=2.0) + gross = uptake(:NO₃, nitrate) + @test [gross, tendency(:P_1, nitrate), tendency(:NH₄, nitrate), tendency(:DOM, nitrate)] ≈ + [1 - exp(-1), 0.95 * gross, 0.0375 * gross, 0.0125 * gross] + + mixed = _fields(; NO₃=10.0, NH₄=10.0, P_1=2.0) + @test uptake(:NO₃, mixed) + uptake(:NH₄, mixed) ≈ + 2 * (1 - exp(-1)) * (10 / 11 * exp(-1) + 10 / 11) + @test [tendency(t, _fields(; Z_1=2.0)) for t in (:Z_1, :NH₄, :DOM)] ≈ [-2.0, 1.0, 1.0] + @test [tendency(t, _fields(; DOM=3.0, H_1=2.0)) for t in (:DOM, :H_1, :NH₄)] ≈ [-3.0, 0.75, 2.25] + + selective_temperature = _controlled(; parameters=(temperature_q10=(P_1=2.0, P_2=1.0),)) + warm = _fields(; NO₃=1.0, T=30.0, P_1=1.0, P_2=1.0) + p1 = selective_temperature(1, 1, 1, _GRID, Val(:P_1), (; time=0.0), warm, aux) + p2 = selective_temperature(1, 1, 1, _GRID, Val(:P_2), (; time=0.0), warm, aux) + @test p1 ≈ 2 * p2 +end diff --git a/test/test_library.jl b/test/test_library.jl index 5b62e426..b1f500a8 100644 --- a/test/test_library.jl +++ b/test/test_library.jl @@ -3,15 +3,20 @@ using Test using ForwardDiff using Agate.Library.Allometry: - consumer_assimilation_matrix_axes, palatability_matrix_allometric_axes, - resolve_diameter_indexed_vector + AllometricParam, SplitPowerLaw, allometric_scaling_power, + palatability_matrix_allometric_axes, + resolve_diameter_indexed_vector, resolve_param using Agate.Library.Nutrients: - frank_tnorm, liebig_minimum, normalized_droop_limitation, quota_uptake_regulation -using Agate.Library.Photosynthesis: geider_light_response, smith_light_limitation + frank_tnorm, inhibited_monod_limitation, liebig_minimum, normalized_droop_limitation, + quota_uptake_regulation +using Agate.Library.Photosynthesis: + exponential_light_limitation, geider_light_response, smith_light_limitation using Agate.Library.Predation: holling_type_ii @testset "Library" begin @test holling_type_ii(1.0, 1.0) == 0.5 + @test exponential_light_limitation(33.0, 33.0) ≈ 1 - exp(-1) + @test inhibited_monod_limitation(1.0, 2.0, 1.0, 0.5) ≈ 0.5 * exp(-1) end @testset "Allometry accepts realized diameter tuples" begin @@ -25,15 +30,20 @@ end consumer_indices=(2,), prey_indices=(1, 2), ) == [1.0 0.5] - @test consumer_assimilation_matrix_axes( - Float64; assimilation_efficiency=[0.2, 0.8], consumer_indices=(2,), prey_indices=(1, 2) - ) == [0.8 0.8] @test_throws ArgumentError resolve_diameter_indexed_vector( Float64, diameters, (true,), 3.0; default=0.0 ) - @test_throws ArgumentError consumer_assimilation_matrix_axes( - Float64; assimilation_efficiency=[0.2, 0.8], consumer_indices=(2,), prey_indices=(3,) - ) +end + +@testset "Split power-law allometry" begin + law = SplitPowerLaw() + c = (prefactor=1.2066, breakpoint=3.0, small_exponent=0.28, large_exponent=-0.15) + at_break = allometric_scaling_power(c.prefactor, c.small_exponent, c.breakpoint) + @test law(c, 1.2) ≈ allometric_scaling_power(c.prefactor, c.small_exponent, 1.2) + @test law(c, c.breakpoint) ≈ at_break + @test law(c, 6.0) ≈ at_break * (6 / c.breakpoint)^(3 * c.large_exponent) + @test resolve_param(Float32, AllometricParam(law; c...), 6.0) isa Float32 + @test_throws ArgumentError law(merge(c, (breakpoint=0.0,)), 6.0) end @testset "Library scalar genericity" begin @@ -41,9 +51,11 @@ end @test Agate.Library.Allometry.allometric_scaling_power(T(1), T(-0.1), T(2)) isa T @test Agate.Library.Nutrients.monod_limitation(T(1), T(0.5)) isa T + @test inhibited_monod_limitation(T(1), T(0.5), T(0.5), T(2)) isa T @test frank_tnorm(T(0.2), T(0.4)) isa T @test frank_tnorm(T(0.2), T(0.4); sharpness=50.0) isa T @test smith_light_limitation(T(50), T(0.1), T(1)) isa T + @test exponential_light_limitation(T(50), T(33)) isa T @test Agate.Library.Mortality.linear_loss(T(1), T(0.1)) isa T @test Agate.Library.Predation.holling_type_ii(T(1), T(0.5)) isa T @test Agate.Library.Remineralization.linear_remineralization(T(1), T(0.1)) isa T diff --git a/test/test_models_construct.jl b/test/test_models_construct.jl index a94c8e22..b4a9354f 100644 --- a/test/test_models_construct.jl +++ b/test/test_models_construct.jl @@ -115,7 +115,7 @@ using Oceananigans.Biogeochemistry: :microzoo_2, ) @test size(named.parameters.palatability_matrix) == (3, 5) - @test size(named.parameters.assimilation_matrix) == (3, 5) + @test named.parameters.assimilation_matrix == fill(Float32(0.32), 3, 5) invalid_size_structures = ( 1, diff --git a/test/test_multistate_process_compilation.jl b/test/test_multistate_process_compilation.jl index b8095088..35c14783 100644 --- a/test/test_multistate_process_compilation.jl +++ b/test/test_multistate_process_compilation.jl @@ -98,6 +98,47 @@ using Agate.Processes: ) end + @testset "fixed-stoichiometry growth products close every growth element" begin + components = ( + P=Plankton(; states=(carbon=:carbon,), reference_state=:carbon), + DIC=Pool(:carbon), DIN=Pool(:nitrogen), DOC=Pool(:carbon), DON=Pool(:nitrogen), + ) + stoichiometry = FixedStoichiometry(; + reference_element=:carbon, bindings=(ratio=(nitrogen=:nitrogen_to_carbon,),) + ) + growth(products) = Growth(; + plankton=:P, reference_resource=:DIC, additional_resources=(nitrogen=:DIN,), + stoichiometry, products, + bindings=(maximum_rate=:maximum_growth_rate, product_fraction=:exudation_fraction), + ) + parameters = ( + maximum_growth_rate=Parameter(0.5), nitrogen_to_carbon=Parameter(0.2), + exudation_fraction=Parameter(0.2), + ) + + products = Products((exudate=(carbon=:DOC, nitrogen=:DON),); stoichiometry) + model = Agate.Construction.construct(ModelDefinition(; + components, processes=(growth=growth(products),), parameters, + )) + test_tendencies( + model, (P=2.0, DIC=10.0, DIN=10.0, DOC=0.0, DON=0.0), + (P=0.8, DIC=-1.0, DIN=-0.2, DOC=0.2, DON=0.04), + ) + + mismatched = FixedStoichiometry(; + reference_element=:carbon, bindings=(ratio=(nitrogen=:other_nitrogen_to_carbon,),) + ) + for invalid in ( + Products(:DOC), + Products((exudate=(carbon=:DOC,),); stoichiometry), + Products((exudate=(carbon=:DOC, nitrogen=:DON),); stoichiometry=mismatched), + ) + @test_throws ArgumentError Agate.Construction.construct(ModelDefinition(; + components, processes=(growth=growth(invalid),), parameters, + )) + end + end + @testset "growth leaves independently prognostic elemental states unchanged" begin definition = ModelDefinition(; components=( diff --git a/test/test_process_compilation.jl b/test/test_process_compilation.jl index 5adbb16d..da6c24c3 100644 --- a/test/test_process_compilation.jl +++ b/test/test_process_compilation.jl @@ -62,3 +62,17 @@ end @test occursin("unrealized targets", message) @test occursin(":not_a_realized_tracer", message) end + +@testset "Growth product routing" begin + components = ( + N=Agate.Components.Pool(:nitrogen), DOM=Agate.Components.Pool(:nitrogen), + NH4=Agate.Components.Pool(:nitrogen), + P=Agate.Components.Plankton(; states=(nitrogen=:nitrogen,), reference_state=:nitrogen, size_structure=[1.0]), + ) + process = Agate.Processes.Growth(; plankton=:P, reference_resource=:N, + bindings=(maximum_rate=:mu, product_fraction=:exudation), + products=Agate.Processes.Products((dissolved=:DOM, inorganic=:NH4); fractions=(inorganic=:share,))) + parameters = (mu=Parameter(2.0), exudation=Parameter(0.2), share=Parameter(0.75)) + bgc = Agate.Construction.construct(ModelDefinition(; components, processes=(growth=process,), parameters)) + test_tendencies(bgc, (N=10.0, DOM=0.0, NH4=0.0, P_1=3.0), (N=-6.0, P_1=4.8, NH4=0.9, DOM=0.3)) +end diff --git a/test/test_processes.jl b/test/test_processes.jl index 301d1082..63cf9562 100644 --- a/test/test_processes.jl +++ b/test/test_processes.jl @@ -3,8 +3,9 @@ using Agate.Components: Plankton, Pool, PlanktonStateRef using Agate.ModelFamilies: default_components, default_processes using Agate.Parameters: ConstantDefault, DerivedDefault, ConstructionParameter, Parameter using Agate.Processes: - AbstractFactor, AbstractFormulation, FactorizedGrowth, Smith, Geider, Monod, - NormalizedDroop, QuotaRegulatedMonod, Liebig, FrankTNorm, Q10, Growth, Light, + AbstractFactor, AbstractFormulation, FactorizedGrowth, Smith, Geider, + ExponentialSaturation, Monod, InhibitedMonod, NormalizedDroop, QuotaRegulatedMonod, + Liebig, FrankTNorm, Q10, Growth, Light, NutrientLimitation, Temperature, NutrientResponse, QuotaResponse, NutrientUptake, FixedStoichiometry, Consumption, Mortality, Products, ModelDefinition, driver_identities, formulation, HeterotrophicConsumption, LinearMortality, @@ -41,6 +42,7 @@ Agate.Processes.factor_value( light = Light(Smith(); driver=:PAR) response = NutrientResponse(Monod(); resource=:N) + inhibited = NutrientResponse(InhibitedMonod(); resource=:N, inhibitor=:A) growth = Growth(; plankton=:P, reference_resource=:N, @@ -48,6 +50,9 @@ Agate.Processes.factor_value( ) @test formulation(light) isa Smith @test formulation(response) isa Monod + @test formulation(inhibited) isa InhibitedMonod + @test formulation(Light(ExponentialSaturation(); driver=:PAR)) isa ExponentialSaturation + @test formulation(Light(Monod(); driver=:PAR)) isa Monod @test participants(growth) == (plankton=(:P,), resource=(:N,)) @@ -182,7 +187,6 @@ Agate.Processes.factor_value( @test_throws ArgumentError canonicalize_model(wrong_element) # Invalid built-in formulation combinations are rejected by their concrete objects or factor contract. - @test_throws MethodError Light(Monod(), :PAR, NamedTuple()) @test_throws MethodError Mortality(Monod(), (:P,), nothing, NamedTuple()) @test_throws ArgumentError canonicalize_model(ModelDefinition(; components=( @@ -195,13 +199,15 @@ end @testset "Built-in parameter domains" begin nodes = ( - FactorizedGrowth(), Smith(), Geider(), Monod(), NormalizedDroop(), - QuotaRegulatedMonod(), FrankTNorm(), Q10(), PreferentialGrazing(), + FactorizedGrowth(), Smith(), Geider(), ExponentialSaturation(), Monod(), + InhibitedMonod(), NormalizedDroop(), QuotaRegulatedMonod(), FrankTNorm(), Q10(:plankton), + PreferentialGrazing(), HeterotrophicConsumption(), LinearMortality(), QuadraticMortality(), LinearRemineralization(), Products((a=:A, b=:B); fractions=(a=:fraction_a,)), FixedStoichiometry(; reference_element=:carbon), ) - expected(node, name) = node isa HeterotrophicConsumption && name === :half_saturation ? :positive : + expected(node, name) = node isa ExponentialSaturation && name === :light_scale ? :positive : + node isa HeterotrophicConsumption && name === :half_saturation ? :positive : name in (:minimum_quota, :maximum_quota, :hill, :sharpness, :q10) ? :positive : name === :reference_temperature ? :finite : name in (:assimilation, :fraction) ? :unit_interval : :nonnegative diff --git a/test/test_recipe_serialization.jl b/test/test_recipe_serialization.jl index 8a25acee..53a066c8 100644 --- a/test/test_recipe_serialization.jl +++ b/test/test_recipe_serialization.jl @@ -1,5 +1,6 @@ using Agate.Construction: decode_recipe, encode_recipe, export_recipe, import_recipe using Agate.ModelFamilies: definition_version +using Agate.Library.Allometry: AllometricParam, SplitPowerLaw using Agate.Models: NiPiZD using OceanBioME: BoxModelGrid using Oceananigans.Biogeochemistry: required_biogeochemical_tracers, biogeochemical_drift_velocity @@ -63,12 +64,13 @@ end "provenance", "content_hash", )) - @test encoded["schema"] == Agate.Construction.recipe_schema() == "agate.model_recipe.v1" + @test encoded["schema"] == Agate.Construction.recipe_schema() == "agate.model_recipe.v0.2" @test encoded["family"] == "NiPiZD" @test encoded["definition_version"] == "0.2.0" @test Set(keys(encoded["realization"])) == Set(( "plankton_pfts", "parameter_overrides", + "setting_overrides", "sinking_tracers", "open_bottom", )) @@ -88,10 +90,20 @@ end @test recipe.parameter_overrides == merge( inputs.parameters, (palatability_matrix=inputs.palatability_matrix,) ) + @test isempty(recipe.setting_overrides) @test !recipe.open_bottom @test recipe.sinking_tracers == inputs.sinking_tracers @test decoded == recipe + split_recipe = Agate.Construction.capture_model_recipe( + family; plankton_pfts=(P=(P=[1.0, 4.0],), Z=(Z=[10.0],)), + parameter_overrides=(maximum_growth_rate=AllometricParam( + SplitPowerLaw(); prefactor=1.2066 / 86400, breakpoint=3.0, + small_exponent=0.28, large_exponent=-0.15, + ),), + ) + @test decode_recipe(encode_recipe(split_recipe)) == split_recipe + mapping_a = (P=(small=[2.0, 1.0], large=[3.0]), Z=(Z=[10.0],)) mapping_b = (Z=(Z=[10.0],), P=(large=[3.0], small=[1.0, 2.0])) overrides_a = (alpha=(small_2=0.3,), maximum_growth_rate=(large_1=1.0e-5, small_1=2.0e-5)) @@ -122,6 +134,7 @@ end microzoo=(:microzoo_1, :microzoo_2), ) @test decoded_manifest == manifest + @test isempty(decoded_manifest.settings) @test decoded_manifest.sinking_tracers.D isa Float32 unsized_recipe = Agate.Construction.ModelRecipe( @@ -129,6 +142,7 @@ end recipe.definition_version, merge(recipe.plankton_pfts, (P=(diat=nothing,),)), recipe.parameter_overrides, + recipe.setting_overrides, recipe.sinking_tracers, recipe.open_bottom, ) @@ -209,6 +223,7 @@ end v"0.2.1", recipe.plankton_pfts, recipe.parameter_overrides, + recipe.setting_overrides, recipe.sinking_tracers, recipe.open_bottom, ) @@ -218,7 +233,7 @@ end (:N, :D, :P_1, :P_2, :Z_1, :Z_2) invalid_schema = modified(encoded) do x - x["schema"] = "agate.model_recipe.invalid" + x["schema"] = "agate.model_recipe.v0.1" end invalid_realization = rehashed(encoded) do x pop!(x["realization"]["plankton_pfts"])