From 5781501bf2b25e8ecf8c023588a4ca9108bc25fb Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Sun, 29 Dec 2024 20:02:44 +0000 Subject: [PATCH 1/6] Work-in-progress updating dimensions and coordinates Simplify and align with Common Data Model Overall design: - dimensions now just have a name and a size - coordinates can optionally be attached to a dimension (currently this are separate Dicts etc) - coordinates are specified by (coordinate) variable name which is a "standard" PALEO variable. Coordinate values then are looked up as needed Changes: - NamedDimension now just contains name, size. Remove FixedCoord (simpler version now in PALEOmodel) - add functions get_dimensions get_dimension set_coordinates! get_coordinates - simplify grids - remove coordinates, which can now be defined as "standard" PALEO variables - remove get_region (similar functionality now in PALEOmodel) - Experimental implementation of Julia CommonDataModel interface in CommonDataModelExt extension. --- Project.toml | 14 +- docs/src/DomainsVariablesFields.md | 1 - ext/CommonDataModelExt/CommonDataModelExt.jl | 93 +++ src/CoordsDims.jl | 203 +------ src/Domain.jl | 92 ++- src/Fields.jl | 229 ++++--- src/Grids.jl | 603 ++++++++++--------- src/Model.jl | 15 +- src/PALEOboxes.jl | 2 +- src/Types.jl | 14 + src/VariableDomain.jl | 18 +- src/reactioncatalog/GridForcings.jl | 2 +- src/reactioncatalog/GridReactions.jl | 58 +- test/runfieldtests.jl | 12 +- 14 files changed, 752 insertions(+), 604 deletions(-) create mode 100644 ext/CommonDataModelExt/CommonDataModelExt.jl diff --git a/Project.toml b/Project.toml index d72dd782..b3742098 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PALEOboxes" uuid = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4" authors = ["Stuart Daines "] -version = "0.21.40" +version = "0.22.0" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" @@ -30,9 +30,16 @@ TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" +[weakdeps] +CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" + +[extensions] +CommonDataModelExt = "CommonDataModel" + [compat] Atomix = "0.1, 1.0" BenchmarkTools = "1.0" +CommonDataModel = "0.3.7" DataFrames = "1.1" DocStringExtensions = "0.8, 0.9" Documenter = "1" @@ -52,12 +59,13 @@ StructArrays = "0.6, 0.7" TestEnv = "1.0" TimerOutputs = "0.5" YAML = "0.4.7" -julia = "1.6" +julia = "1.10" [extras] +CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Documenter", "Logging", "Test"] +test = ["CommonDataModel", "Documenter", "Logging", "Test"] diff --git a/docs/src/DomainsVariablesFields.md b/docs/src/DomainsVariablesFields.md index 91c479c0..72503bba 100644 --- a/docs/src/DomainsVariablesFields.md +++ b/docs/src/DomainsVariablesFields.md @@ -148,7 +148,6 @@ Examples: ```@docs Field get_field -wrap_field ``` ## Spaces diff --git a/ext/CommonDataModelExt/CommonDataModelExt.jl b/ext/CommonDataModelExt/CommonDataModelExt.jl new file mode 100644 index 00000000..c2bdcfc4 --- /dev/null +++ b/ext/CommonDataModelExt/CommonDataModelExt.jl @@ -0,0 +1,93 @@ +module CommonDataModelExt + +import PALEOboxes as PB +import CommonDataModel as CDM + +##################################### +# Wrapper types +###################################### + +struct ModelCDM <: CDM.AbstractDataset + model::PB.Model + modeldata::Union{Nothing, PB.AbstractModelData} +end + +PB.CDModel(model::PB.Model) = ModelCDM(model, nothing) +PB.CDModel(modeldata::PB.AbstractModelData) = ModelCDM(modeldata.model, modeldata) + +struct DomainCDM <: CDM.AbstractDataset + domain::PB.Domain + modeldata::Union{Nothing, PB.AbstractModelData} +end + +PB.CDModel(domain::PB.Domain) = DomainCDM(domain, nothing) + + +struct VariableDomainCDM <: CDM.AbstractDataset + variabledomain::PB.VariableDomain + modeldata::Union{Nothing, PB.AbstractModelData} + data +end + +PB.CDModel(variabledomain::PB.VariableDomain) = VariableDomainCDM(variabledomain, nothing, nothing) + +###################################### +# Model +###################################### + +# iterable with all group names +CDM.groupnames(m::ModelCDM) = [d.name for d in m.model.domains] + +CDM.group(m::ModelCDM, name::AbstractString) = DomainCDM(PB.get_domain(m.model, name; allow_not_found=false), m.modeldata) + +############################################### +# Domain +################################################ + +CDM.name(d::DomainCDM) = d.domain.name + +# TODO +# parentdataset(d::DomainCDM) + +# returns a list of variable names as strings +Base.keys(d::DomainCDM) = [v.name for v in PB.get_variables(d.domain)] + +function CDM.variable(d::DomainCDM, varname::AbstractString) + variabledomain = PB.get_variable(d.domain, varname; allow_not_found=false) + if isnothing(d.modeldata) + data = nothing + else + data = PB.get_data(variabledomain, d.modeldata) + end + return VariableDomainCDM(variabledomain, d.modeldata, data) +end + +CDM.dimnames(d::DomainCDM) = [nd.name for nd in PB.get_dimensions(d.domain)] + +CDM.dim(d::DomainCDM, name::AbstractString) = PB.get_dimension(d.domain, name).size + + +############################################### +# VariableDomain +################################################ + +CDM.name(v::VariableDomainCDM) = v.variabledomain.name + +CDM.dataset(v::VariableDomainCDM) = DomainCDM(v.variabledomain.domain, v.modeldata) + +CDM.dimnames(v::VariableDomainCDM) = [nd.name for nd in PB.get_dimensions(v.variabledomain)] + +Base.ndims(v::VariableDomainCDM) = length(CDM.dimnames(v)) + +Base.size(v::VariableDomainCDM) = (nd.size for nd in PB.get_dimensions(v.variabledomain)) + +CDM.attribnames(v::VariableDomainCDM) = keys(v.variabledomain.attributes) + +CDM.attrib(v::VariableDomainCDM, name::Symbol) = v.variabledomain.attributes[name] + +Base.getindex(v::VariableDomainCDM, indices...) = Base.getindex(v.data, indices...) + +Base.eltype(v::VariableDomainCDM) = Base.eltype(v.data) + + +end # module \ No newline at end of file diff --git a/src/CoordsDims.jl b/src/CoordsDims.jl index cbe9d930..80bb662b 100644 --- a/src/CoordsDims.jl +++ b/src/CoordsDims.jl @@ -1,199 +1,56 @@ -################################ -# Coordinates -################################# -""" - FixedCoord -A fixed (state independent) coordinate -""" -mutable struct FixedCoord - name::String - values::Vector{Float64} - attributes::Dict{Symbol, Any} -end +################################################# +# Dimensions +##################################################### """ - append_units(name::AbstractString, attributes) -> "name (units)" + NamedDimension(name, size) -Utility function to append variable units string to a variable name for display. +A named dimension """ -function append_units(name::AbstractString, attributes::Dict{Symbol, Any}) - units = get(attributes, :units, "") - if isempty(units) - return name - else - return name*" ($units)" - end +struct NamedDimension + name::String + size::Int64 end -append_units(name::AbstractString, attributes::Nothing) = name +function Base.show(io::IO, nd::NamedDimension) + print(io, "NamedDimension(name=", nd.name, ", size=", nd.size, ")") + return nothing +end """ - build_coords_edges(coords_vec::Vector{FixedCoord}) -> Vector{Float64} - -Build a vector of coordinate edges (length `n+1``) from `coords_vec`, assuming the PALEO -convention that `coords_vec` contains three elements with -cell midpoints, lower edges, upper edges each of length `n`, in that order. + function get_dimensions(obj) -> Vector{NamedDimension} -Falls back to just returning the first entry in `coords_vec` for other cases. +Get all dimensions for PALEO object `obj` """ -function build_coords_edges(coords_vec::Vector{FixedCoord}) - - if length(coords_vec) == 1 || length(coords_vec) > 3 - # 1 coordinate or something we don't understand - take first - co = first(coords_vec) - co_values = co.values - co_label = append_units(co.name, co.attributes) - elseif length(coords_vec) in (2, 3) - # 2 coordinates assume lower, upper edges - # 3 coordinates assume mid, lower, upper - co_lower = coords_vec[end-1] - co_upper = coords_vec[end] - co_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes) - first(co_lower.values) < first(co_upper.values) || - @warn "build_coords_edges: $co_label co_lower is > co_upper - check model grid" - if co_lower.values[end] > co_lower.values[1] # ascending order - co_lower.values[2:end] == co_upper.values[1:end-1] || - @warn "build_coords_edges: $co_label lower and upper edges don't match" - co_values = [co_lower.values; co_upper.values[end]] - else # descending order - co_lower.values[1:end-1] == co_upper.values[2:end] || - @warn "build_coords_edges: $co_label lower and upper edges don't match" - co_values = [co_upper.values[1]; co_lower.values] - end - - end - - return co_values, co_label -end +function get_dimensions end -"guess coordinate edges from midpoints, assuming uniform spacing" -function guess_coords_edges(x_midpoints) - first_x = x_midpoints[1] - 0.5*(x_midpoints[2] - x_midpoints[1]) - last_x = x_midpoints[end] + 0.5*(x_midpoints[end] - x_midpoints[end-1]) - return [first_x; 0.5.*(x_midpoints[1:end-1] .+ x_midpoints[2:end]); last_x] -end - - -function get_region(fc::FixedCoord, indices::AbstractVector) - return FixedCoord(fc.name, fc.values[indices], fc.attributes) -end - -function get_region(fcv::Vector{FixedCoord}, indices::AbstractVector) - return [FixedCoord(fc.name, fc.values[indices], fc.attributes) for fc in fcv] -end - - -"find indices of coord from first before range[1] to first after range[2]" -function find_indices(coord::AbstractVector, range) - length(range) == 2 || - throw(ArgumentError("find_indices: length(range) != 2 $range")) - - idxstart = findlast(t -> t<=range[1], coord) - isnothing(idxstart) && (idxstart = 1) - - idxend = findfirst(t -> t>=range[2], coord) - isnothing(idxend) && (idxend = length(coord)) - - return idxstart:idxend, (coord[idxstart], coord[idxend]) -end - -"find indices of coord nearest val" -function find_indices(coord::AbstractVector, val::Real) - idx = 1 - for i in 1:length(coord) - if abs(coord[i] - val) < abs(coord[idx] - val) - idx = i - end - end - - return [idx], coord[idx] -end +""" + function get_dimension(obj, dimname) -> NamedDimension -################################################# -# Dimensions -##################################################### +Get all dimension `dimname` for PALEO object `obj` +""" +function get_dimension end """ - NamedDimension + function set_coordinates!(obj, dimname, coordinates::Vector{String}) -A named dimension, with optional attached fixed coordinates `coords` +Set coordinates attached to `dimname` for PALEO object `obj` -PALEO convention is that where possible `coords` contains three elements, for cell +PALEO convention is that where possible `coordinates` contains three elements, for cell midpoints, lower edges, upper edges, in that order. """ -mutable struct NamedDimension - name::String - size::Int64 - coords::Vector{FixedCoord} # may be empty -end - -"create from size only (no coords)" -function NamedDimension(name, size::Integer) - return NamedDimension( - name, - size, - FixedCoord[], - ) -end -"create from coord mid-points" -function NamedDimension(name, coord::AbstractVector) - return NamedDimension( - name, - length(coord), - [ - FixedCoord(name, coord, Dict{Symbol, Any}()), - ] - ) -end - -"create from coord mid-points and edges" -function NamedDimension(name, coord::AbstractVector, coord_edges::AbstractVector) - if coord[end] > coord[1] - # ascending order - coord_lower = coord_edges[1:end-1] - coord_upper = coord_edges[2:end] - else - # descending order - coord_lower = coord_edges[2:end] - coord_upper = coord_edges[1:end-1] - end - return NamedDimension( - name, - length(coord), - [ - FixedCoord(name, coord, Dict{Symbol, Any}()), - FixedCoord(name*"_lower", coord_lower, Dict{Symbol, Any}()), - FixedCoord(name*"_upper", coord_upper, Dict{Symbol, Any}()), - ] - ) -end - -function get_region(nd::NamedDimension, indices::AbstractVector) - return NamedDimension(nd.name, length(indices), get_region(nd.coords, indices)) -end +function set_coordinates! end """ - build_coords_edges(nd::NamedDimension) -> Vector{Float64} + function get_coordinates(obj, dimname) -> coordinates::Vector{String} -Call [`build_coords_edges`](@ref)(nd.coords), or fallback to just returning indices -if no coords present. -""" -function build_coords_edges(nd::NamedDimension) - if !isempty(nd.coords) - return build_coords_edges(nd.coords) - else - @warn "no coords for NamedDimension $(nd.name), returning indices" - return collect(1:nd.size), nd.name*" (indices)" - end -end +Get coordinates (if any) attached to `dimname` for PALEO object `obj` -function Base.show(io::IO, nd::NamedDimension) - print(io, "NamedDimension(name=", nd.name, ", size=", nd.size, ", coords=(") - join(io, [c.name for c in nd.coords], ", ") - print(io, "))") - return nothing -end +PALEO convention is that where possible `coordinates` contains three elements, for cell +midpoints, lower edges, upper edges, in that order. +""" +function get_coordinates end diff --git a/src/Domain.jl b/src/Domain.jl index 1c9e13d1..6f33fe25 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -14,23 +14,23 @@ Base.@kwdef mutable struct Domain <: AbstractDomain name::String ID::Int data_dims::Vector{NamedDimension} = Vector{NamedDimension}() + data_dims_coordinates::Dict{String, Vector{String}} = Dict{String, Vector{String}}() parameters::Dict{String, Any} - grid::Union{Nothing, AbstractMesh} = nothing - + grid::Union{Nothing, AbstractMesh} = nothing reactions::Vector{AbstractReaction} = Vector{AbstractReaction}() variables::Dict{String, VariableDomain} = Dict{String, VariableDomain}() end """ - set_data_dimension!(domain::Domain, dim::NamedDimension; allow_exists=false) + set_data_dimension!(domain::Domain, dim::NamedDimension [, coordinates], ; allow_exists=false) -Define a Domain data dimension as a [`NamedDimension`](@ref) +Define a Domain data dimension as a [`NamedDimension`](@ref), optionally attaching a Vector of coordinate names. Variables may then specify data dimensions as a list of names using the `:data_dims` Variable Attribute. """ -function set_data_dimension!(domain::Domain, dim::NamedDimension; allow_exists=false) +function set_data_dimension!(domain::Domain, dim::NamedDimension, coordinates::Vector{String} = String[]; allow_exists=false) - @info "set_data_dimension!: setting Domain '$(domain.name)' data dimension '$dim'" + @info "set_data_dimension!: setting Domain '$(domain.name)' data dimension '$dim' coordinates $coordinates" idx = findfirst(d -> d.name==dim.name, domain.data_dims) @@ -44,19 +44,42 @@ function set_data_dimension!(domain::Domain, dim::NamedDimension; allow_exists=f domain.data_dims[idx] = dim end + delete!(domain.data_dims_coordinates, dim.name) + if !isempty(coordinates) + domain.data_dims_coordinates[dim.name] = coordinates + end + return nothing end -has_data_dimension(domain::Domain, dimname::AbstractString) = - !isnothing(findfirst(d -> d.name==dimname, domain.data_dims)) - function get_data_dimension(domain::Domain, dimname::AbstractString) idx = findfirst(d -> d.name==dimname, domain.data_dims) !isnothing(idx) || - error("Domain $(domain.name) has no dimension='$dimname' (available dimensions: $(domain.data_dims)") + error("Domain $(domain.name) has no data_dimension='$dimname' (available data_dimensions: $(domain.data_dims)") return domain.data_dims[idx] end + +function get_dimensions(domain::Domain) + spatial_dims = isnothing(domain.grid) ? NamedDimension[] : get_dimensions(domain.grid) + return vcat(spatial_dims, domain.data_dims) +end + +function get_dimension(domain::Domain, dimname::AbstractString) + ad = get_dimensions(domain) + idx = findfirst(d -> d.name==dimname, ad) + !isnothing(idx) || + error("Domain $(domain.name) has no dimension='$dimname' (available dimensions: $ad") + return ad[idx] +end + +function set_coordinates!(domain::Domain, dimname, coordinates::Vector{String}) + idx = findfirst(nd -> nd.name == dimname, domain.data_dims) + !isnothing(idx) || error("set_coordinates! domain $(domain.name) has no data dimension $dimname") + domain.data_dims_coordinates[dimname] = coordinates + return nothing +end + function get_length(domain::Domain) if isnothing(domain.grid) return 1 # scalar Domain @@ -65,6 +88,7 @@ function get_length(domain::Domain) end end + "Get number of Domain variables" function get_num_variables(domain::Domain) return length(domain.variables) @@ -536,28 +560,46 @@ end # Pretty printing ############################' -"compact form" + # compact form function Base.show(io::IO, domain::Domain) print(io, "Domain(name='", domain.name, "')") end -"multiline form" -function Base.show(io::IO, ::MIME"text/plain", domain::Domain) +# multiline form +function Base.show(io::IO, ::MIME"text/plain", domain::Domain; show_reactions=true, show_variables=true) println(io, "Domain") println(io, " name='$(domain.name)'") println(io, " ID=$(domain.ID)") - println(io, " data_dims=", domain.data_dims) - println(io, " grid=", isnothing(domain.grid) ? "" : domain.grid) - println(io, " reactions:") - for r in domain.reactions - println(io, " ", r) + println(io, " data_dims:") + for nd in domain.data_dims + println(io, " ", nd) + if haskey(domain.data_dims_coordinates, nd.name) + println(io, " coordinates: ", domain.data_dims_coordinates[nd.name]) + end end - println(io, " variables (VariableDomPropDep):") - for var in sort(get_variables(domain, v -> v isa VariableDomPropDep), by = v -> v.name) - println(io, " ", var) + println(io, " grid:") + if !isnothing(domain.grid) + iogrid = IOBuffer() + show(iogrid, MIME"text/plain"(), domain.grid) + seekstart(iogrid) + for line in eachline(iogrid) + println(io, " ", line) + end + end + if show_reactions + println(io, " reactions:") + for r in domain.reactions + println(io, " ", r) + end end - println(io, " variables (VariableDomContribTarget):") - for var in sort(get_variables(domain, v -> v isa VariableDomContribTarget), by = v -> v.name) - println(io, " ", var) + if show_variables + println(io, " variables (VariableDomPropDep):") + for var in sort(get_variables(domain, v -> v isa VariableDomPropDep), by = v -> v.name) + println(io, " ", var) + end + println(io, " variables (VariableDomContribTarget):") + for var in sort(get_variables(domain, v -> v isa VariableDomContribTarget), by = v -> v.name) + println(io, " ", var) + end end end @@ -575,7 +617,7 @@ Show table of Domain Variables. Optionally get variable links, data. """ function show_variables( domain::Domain; - attributes=[:units, :vfunction, :space, :field_data, :description], + attributes=[:units, :vfunction, :space, :data_dims, :field_data, :description], filter = attrb->true, showlinks=false, modeldata=nothing diff --git a/src/Fields.jl b/src/Fields.jl index 27c69b6e..2690e5e6 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -6,7 +6,7 @@ """ AbstractSpace -Defines a function space within a Domain, on a mesh defined by a Grid +Defines a function Space within a Domain, on a mesh defined by a Grid """ AbstractSpace @@ -59,7 +59,15 @@ end ################################################################# """ - internal_size(::Type{<:AbstractSpace}, mesh::AbstractMesh; [subdomain=""] [space=:cell]) -> NTuple{ndims, Int} + has_internal_cartesian(mesh::AbstractMesh, Space::Type{<:AbstractSpace})::Bool + +`true` if `Mesh` uses different internal and external (cartesian) spatial array representations for `Space` +""" +function has_internal_cartesian end + + +""" + internal_size(Space::Type{<:AbstractSpace}, mesh::AbstractMesh; [subdomain=""] [space=:cell]) -> NTuple{ndims, Int} Array size to use for model Variables. @@ -71,7 +79,7 @@ All `AbstractMesh` concrete subtypes (UnstructuredVectorGrid, CartesianLinearGri function internal_size end """ - cartesian_size(mesh::AbstractMesh) -> NTuple{ndims, Int} + cartesian_size(Space::Type{<:AbstractSpace}, mesh::AbstractMesh) -> NTuple{ndims, Int} Optional (only regular Cartesian grids should implement this method): Array size of Cartesian Domain. @@ -80,12 +88,13 @@ eg to a Vector for internal model Variables. """ function cartesian_size end + """ spatial_size(::Type{<:AbstractSpace}, mesh::AbstractMesh) -> NTuple{ndims, Int} Array size for given Space and mesh. """ -spatial_size(space::Type{<:AbstractSpace}, mesh) = internal_size(space, mesh) +spatial_size(Space::Type{<:AbstractSpace}, mesh) = internal_size(Space, mesh) ################################################################ @@ -111,24 +120,24 @@ If the subtype has a representation as components, it should implement: [`num_components`](@ref), [`get_components`](@ref) If the subtype needs to provide a thread-safe atomic addition operation eg to provide scalar accumulator variables for Domain totals with a tiled model, -it should implement [`atomic_add!`](@ref) for the field `data_type`. +it should implement [`atomic_add!`](@ref) for the field values. """ AbstractData """ allocate_values( - field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, space::Type{<:AbstractSpace}, spatial_size::Tuple; + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, Space::Type{<:AbstractSpace}, spatial_size::Tuple; allocatenans::Bool, - ) -> values + ) -> values::V -allocate `Field.values` (eg an Array) for `field_data` with dimensions defined by `spatial_size` and `data_dims` +allocate `Field.values::V` (eg an Array of elements with Type `data_type`) for `FieldData` with dimensions defined by `spatial_size` and `data_dims` """ function allocate_values( - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, spatial_size::Tuple{Integer, Vararg{Integer}}, # an NTuple with at least one element allocatenans::Bool, ) @@ -137,17 +146,17 @@ end """ check_values( existing_values, - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, spatial_size::Tuple{Integer, Vararg{Integer}} ) Check `existing_values` is of suitable type, size etc for use as `Field.values`, throw exception if not. """ function check_values( - existing_values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, space::Type{<:AbstractSpace}, spatial_size::Tuple{Integer, Vararg{Integer}}, + existing_values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, data_type, Space::Type{<:AbstractSpace}, spatial_size::Tuple{Integer, Vararg{Integer}}, ) end @@ -168,14 +177,14 @@ end """ init_values!( - values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, + values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, init_value::Symbol, attribv::VariableBase, convertfn, convertvalues, cellrange, info::NTuple{3, String} ) Initialize `values` at model start to `init_value` over region `cellrange` using information from Variable `attribv` attributes, scaled by `convertfn` and `convertvalues`. -Optional: only required if this `field_data` type is used for a model (state) Variable that requires initialisation. +Optional: only required if this `FieldData` type is used for a model (state) Variable that requires initialisation. Arguments: - `values`: data to be zeroed @@ -189,28 +198,28 @@ Arguments: """ function init_values!( - values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, + values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, init_value::Symbol, attribv::VariableBase, convertfn, convertvalues, cellrange, info::NTuple{3, String} ) end """ - zero_values!(values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, cellrange) + zero_values!(values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, cellrange) Set `values` over spatial region `cellrange` to zero at start of main loop """ -function zero_values!(values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, cellrange) +function zero_values!(values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, cellrange) end """ - field_single_element(field_data::Type{<:AbstractData}, N) -> Bool + field_single_element(FieldData::Type{<:AbstractData}, N) -> Bool -Return true if `field_data` with length(data_dims) = N is represented with +Return true if `FieldData` with length(data_dims) = N is represented with a single value that can be accessed as [] (used to optimize FieldRecord storage). -Default is probably OK, unless `field_data` uses a Vector of values per element. +Default is probably OK, unless `FieldData` uses a Vector of values per element. """ -function field_single_element(field_data::Type{<:AbstractData}, N) +function field_single_element(FieldData::Type{<:AbstractData}, N) if N == 0 return true else @@ -220,16 +229,16 @@ end """ dof_values( - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, mesh, cellrange ) -> dof::Int -Return degrees-of-freedom for `field_data` over spatial region `cellrange`. +Return degrees-of-freedom for `FieldData` over spatial region `cellrange`. """ -function dof_values(field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, mesh, cellrange) +function dof_values(FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, mesh, cellrange) end """ @@ -237,9 +246,9 @@ end dest, doff, values, - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, cellrange ) -> num_copied::Int @@ -247,17 +256,17 @@ Copy Field.values `values` from spatial region defined by `cellrange`, to `dest` Number of values over whole Domain should equal degrees-of-freedom returned by [`dof_values`](@ref) -Required if this `field_data` type needs to provide values for a numerical solver. +Required if this `FieldData` type needs to provide values for a numerical solver. """ -function copyfieldto!(dest, doff, values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, cellrange) +function copyfieldto!(dest, doff, values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, cellrange) end """ copytofield!( values, - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, cellrange, src, soff @@ -267,25 +276,25 @@ Copy from `src` Array starting at index `soff` to Field.values `values` for spat Number of values over whole Domain should equal degrees-of-freedom returned by [`dof_values`](@ref) -Required if this `field_data` type needs to provide values for a numerical solver. +Required if this `FieldData` type needs to provide values for a numerical solver. """ -function copytofield!(values, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, cellrange, src, soff) +function copytofield!(values, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, cellrange, src, soff) end """ - add_field!(dest, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, a, cellrange, src) + add_field!(dest, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, a, cellrange, src) Implement `dest += a*src` where `dest`, `src` are Field.values, `a` is a number, over region defined by `cellrange` """ -function add_field!(dest, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, a, cellrange, src) +function add_field!(dest, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, a, cellrange, src) end """ add_field_vec!( dest, - field_data::Type{<:AbstractData}, + FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, - space::Type{<:AbstractSpace}, + Space::Type{<:AbstractSpace}, a, cellrange, src, @@ -299,27 +308,27 @@ Returns number of elements of `src` used. See [`copytofield!`](@ref), [`copyfieldto!`](@ref) for the relationship between Array `src` and Field values `dest`. """ -function add_field_vec!(dest, field_data::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space::Type{<:AbstractSpace}, a, cellrange, src, soff) +function add_field_vec!(dest, FieldData::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space::Type{<:AbstractSpace}, a, cellrange, src, soff) end """ - num_components(field_data::Type{<:AbstractData}) -> Int + num_components(FieldData::Type{<:AbstractData}) -> Int -get number of components (optional - implement if `field_data` has a representation as components) +get number of components (optional - implement if `FieldData` has a representation as components) """ -function num_components(field_data::Type{<:AbstractData}) end +function num_components(FieldData::Type{<:AbstractData}) end """ - get_components(values, field_data::Type{<:AbstractData}) -> Vector + get_components(values, FieldData::Type{<:AbstractData}) -> Vector Convert Field `values` to a Vector of components -(optional - implement if `field_data` has a representation as components) +(optional - implement if `FieldData` has a representation as components) """ -function get_components(values, field_data::Type{<:AbstractData}) end +function get_components(values, FieldData::Type{<:AbstractData}) end "Optional: sanitize `values` for storing as model output. Default implementation is usually OK - only implement for custom types that should be converted to standard types for storage" -get_values_output(values, data_type::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, space, mesh) = values +get_values_output(values, data_type::Type{<:AbstractData}, data_dims::Tuple{Vararg{NamedDimension}}, Space, mesh) = values @@ -343,19 +352,68 @@ end """ - Field{D <: AbstractData, S <: AbstractSpace, V, N, M} + Field{FieldData <: AbstractData, Space <: AbstractSpace, V, N, Mesh <: AbstractMeshOrNothing} + Field(FieldData, data_dims::NTuple{N, NamedDimensions}, data_type, Space, mesh::AbstractMeshOrNothing; allocatenans) + Field(existing_values, FieldData, data_dims::NTuple{N, NamedDimension}, data_type::Union{DataType, Missing}, Space, mesh::AbstractMeshOrNothing) -A Field of `values::V` of data type `D` defined on function space `S` over `mesh::M` +A Field of `values::V` of type `FieldData <: AbstractData` defined on function space `Space` over grid of type `Mesh` and (optionally) with `N` `data_dims::NTuple{N, NamedDimensions}`. + +`values::V` is usually Array or similar of elements with Type `data_type` and dimensions defined by `Space`, `Mesh` and `data_dims`, +and is either created by [`allocate_values`](@ref) or supplied by `existing_values`. """ -struct Field{D <: AbstractData, S <: AbstractSpace, V, N, M} +struct Field{FieldData <: AbstractData, Space <: AbstractSpace, V, N, Mesh <: AbstractMeshOrNothing} values::V data_dims::NTuple{N, NamedDimension} - mesh::M + mesh::Mesh end -field_data(field::Field{D, S, V, N, M}) where {D, S, V, N, M} = D -space(field::Field{D, S, V, N, M}) where {D, S, V, N, M} = S +# create a new Field, allocating `values` data arrays +function Field( + FieldData::Type, data_dims::NTuple{N, NamedDimension}, data_type::Type, Space::Type{<:AbstractSpace}, mesh::AbstractMeshOrNothing; + allocatenans +) where {N} + v = allocate_values( + FieldData, data_dims, data_type, Space, spatial_size(Space, mesh); + allocatenans, + ) + + return Field{FieldData, Space, typeof(v), N, typeof(mesh)}(v, data_dims, mesh) +end + +# create a new Field, containing supplied `existing_values` data arrays +function Field( + existing_values, FieldData::Type, data_dims::NTuple{N, NamedDimension}, data_type::Union{DataType, Missing}, Space::Type{<:AbstractSpace}, mesh::AbstractMeshOrNothing +) where {N} + check_values( + existing_values, FieldData, data_dims, data_type, Space, spatial_size(Space, mesh), + ) + return Field{FieldData, Space, typeof(existing_values), N, typeof(mesh)}(existing_values, data_dims, mesh) +end + +field_data(field::Field{FieldData, Space, V, N, Mesh}) where {FieldData, Space, V, N, Mesh} = FieldData +space(field::Field{FieldData, Space, V, N, Mesh}) where {FieldData, Space, V, N, Mesh} = Space + +function get_dimensions(f::Field; expand_cartesian=false) + dims = NamedDimension[nd for nd in get_dimensions(f.mesh, space(f); expand_cartesian)] + append!(dims, f.data_dims) + return dims +end + +# use default compact form Base.show + +# multiline form +function Base.show(io::IO, ::MIME"text/plain", f::Field{FieldData, Space, V, N, Mesh}) where {FieldData, Space, V, N, Mesh} + println(io, "Field{", FieldData, ", ", Space, ", ", V, ", ", N, ", ", Mesh, "}") + println(io, " values: ", f.values) + println(io, " data_dims: ", f.data_dims) + println(io, " mesh: ", f.mesh) + println(io, " dimensions:") + for nd in get_dimensions(f) + println(io, " ", nd) + end + return nothing +end """ get_field(obj, ...) -> Field @@ -371,100 +429,79 @@ Add Field or Field to PALEO object `obj` """ function add_field! end -"create a new Field, allocating `values` data arrays" -function allocate_field( - field_data::Type, data_dims::NTuple{N, NamedDimension}, data_type::Type, space::Type{<:AbstractSpace}, mesh; - allocatenans -) where {N} - v = allocate_values( - field_data, data_dims, data_type, space, spatial_size(space, mesh); - allocatenans, - ) - - return Field{field_data, space, typeof(v), N, typeof(mesh)}(v, data_dims, mesh) -end -"create a new Field, containing supplied `existing_values` data arrays" -function wrap_field( - existing_values, field_data::Type, data_dims::NTuple{N, NamedDimension}, data_type::Union{DataType, Missing}, space::Type{<:AbstractSpace}, mesh -) where {N} - check_values( - existing_values, field_data, data_dims, data_type, space, spatial_size(space, mesh), - ) - return Field{field_data, space, typeof(existing_values), N, typeof(mesh)}(existing_values, data_dims, mesh) -end "zero out `field::Field` over region defined by `cellrange`" -function zero_field!(field::Field{D, S, V, N, M}, cellrange) where {D, S, V, N, M} - zero_values!(field.values, D, field.data_dims, S, cellrange) +function zero_field!(field::Field{FieldData, Space, V, N, Mesh}, cellrange) where {FieldData, Space, V, N, Mesh} + zero_values!(field.values, FieldData, field.data_dims, Space, cellrange) end "initialize `field::Field` to `init_value` (`:initial_value` or `:norm_value`) from Variable `attribv` attributes, over region defined by `cellrange`. Optionally calculate transformed initial values from supplied `convertfn` and `convertvalues`" function init_field!( - field::Field{D, S, V, N, M}, init_value::Symbol, attribv::VariableBase, convertfn, convertvalues, cellrange, info, -) where {D, S, V, N, M} - init_values!(field.values, D, field.data_dims, S, init_value, attribv, convertfn, convertvalues, cellrange, info) + field::Field{FieldData, Space, V, N, Mesh}, init_value::Symbol, attribv::VariableBase, convertfn, convertvalues, cellrange, info, +) where {FieldData, Space, V, N, Mesh} + init_values!(field.values, FieldData, field.data_dims, Space, init_value, attribv, convertfn, convertvalues, cellrange, info) end "calculate number of degrees-of-freedom for `field::Field` over region defined by `cellrange`" -function dof_field(field::Field{D, S, V, N, M}, cellrange) where {D, S, V, N, M} - return dof_values(D, field.data_dims, S, field.mesh, cellrange) +function dof_field(field::Field{FieldData, Space, V, N, Mesh}, cellrange) where {FieldData, Space, V, N, Mesh} + return dof_values(FieldData, field.data_dims, Space, field.mesh, cellrange) end "copy `src::Field`` to `dest::Vector`, optionally restricting to region defined by `cellrange`" -function Base.copyto!(dest, doff, src::Field{D, S, V, N, M}, cellrange) where {D, S, V, N, M} - return copyfieldto!(dest, doff, src.values, D, src.data_dims, S, cellrange) +function Base.copyto!(dest, doff, src::Field{FieldData, Space, V, N, Mesh}, cellrange) where {FieldData, Space, V, N, Mesh} + return copyfieldto!(dest, doff, src.values, FieldData, src.data_dims, Space, cellrange) end "copy `src::Vector`` to `dest::Field`, optionally restricting to region defined by `cellrange`" -function Base.copyto!(dest::Field{D, S, V, N, M}, cellrange, src, soff) where {D, S, V, N, M} - return copytofield!(dest.values, D, dest.data_dims, S, cellrange, src, soff) +function Base.copyto!(dest::Field{FieldData, Space, V, N, Mesh}, cellrange, src, soff) where {FieldData, Space, V, N, Mesh} + return copytofield!(dest.values, FieldData, dest.data_dims, Space, cellrange, src, soff) end "Calculate `dest::Field = dest::Field + a * src::Field`, optionally restricting to region defined by `cellrange`" -function add_field!(dest::Field{D, S, V, N, M}, a, cellrange, src::Field{D, S, V, N, M}) where {D, S, V, N, M} - return add_field!(dest.values, D, dest.data_dims, S, a, cellrange, src.values) +function add_field!(dest::Field{FieldData, Space, V, N, Mesh}, a, cellrange, src::Field{FieldData, Space, V, N, Mesh}) where {FieldData, Space, V, N, Mesh} + return add_field!(dest.values, FieldData, dest.data_dims, Space, a, cellrange, src.values) end -function add_field_vec!(dest::Field{D, S, V, N, M}, a, cellrange, srcvalues::AbstractVector, soff) where {D, S, V, N, M} - return add_field_vec!(dest.values, D, dest.data_dims, S, a, cellrange, srcvalues, soff) +function add_field_vec!(dest::Field{FieldData, Space, V, N, Mesh}, a, cellrange, srcvalues::AbstractVector, soff) where {FieldData, Space, V, N, Mesh} + return add_field_vec!(dest.values, FieldData, dest.data_dims, Space, a, cellrange, srcvalues, soff) end "sanitized version of `values`, suitable for storing as output" -function get_values_output(field::Field{D, S, V, N, M}) where {D, S, V, N, M} - return get_values_output(field.values, D, field.data_dims, S, field.mesh) +function get_values_output(field::Field{FieldData, Space, V, N, Mesh}) where {FieldData, Space, V, N, Mesh} + return get_values_output(field.values, FieldData, field.data_dims, Space, field.mesh) end # get values from `linkvar_field`, optionally applying view defined by `linksubdomain` function create_accessor( - output_data::Union{Type{D}, Type{UndefinedData}}, + output_data::Union{Type{FieldData}, Type{UndefinedData}}, # TODO - check space, data_dims, - linkvar_field::Field{D, S, V, N, M}, linksubdomain::Union{Nothing, AbstractSubdomain}; + linkvar_field::Field{FieldData, Space, V, N, Mesh}, linksubdomain::Union{Nothing, AbstractSubdomain}; forceview, components, -) where {D, S, V, N, M} +) where {FieldData, Space, V, N, Mesh} # create accessor if isnothing(linksubdomain) if forceview if components - var_accessor = [view(vc, 1:length(vc)) for vc in get_components(linkvar_field.values, D)] + var_accessor = [view(vc, 1:length(vc)) for vc in get_components(linkvar_field.values, FieldData)] else var_accessor = view(linkvar_field.values, 1:length(linkvar_field.values)) end else if components - var_accessor = get_components(linkvar_field.values, D) + var_accessor = get_components(linkvar_field.values, FieldData) else var_accessor = linkvar_field.values end end else if components - var_accessor = [Grids.subdomain_view(vc, linksubdomain) for vc in get_components(linkvar_field.values, D)] + var_accessor = [Grids.subdomain_view(vc, linksubdomain) for vc in get_components(linkvar_field.values, FieldData)] else var_accessor = Grids.subdomain_view(linkvar_field.values, linksubdomain) end diff --git a/src/Grids.jl b/src/Grids.jl index 26f00bf3..11f61d5f 100644 --- a/src/Grids.jl +++ b/src/Grids.jl @@ -129,11 +129,23 @@ Defines additional geometric and topological information for [`PB.Domain`](@ref) Concrete subtypes should implement methods: +[`available_spaces`](@ref) tuple of [`PB.AbstractSpace`](@ref) types supported. + +[`PB.has_internal_cartesian`](@ref) `true` if subtype uses an internal spatial representation that differs from +the external (cartesian) representation. + +[`PB.get_dimensions`](@ref) + [`PB.internal_size`](@ref), optionally [`PB.cartesian_size`](@ref) [`PB.Grids.set_subdomain!`](@ref), [`PB.Grids.get_subdomain`](@ref) -[`PB.Grids.create_default_cellrange`](@ref), [`PB.Grids.get_region`](@ref) +[`PB.Grids.create_default_cellrange`](@ref) + +# Internal and cartesian representations +If a subtype uses an internal representation that differs from the external (cartesian) representation, +it should define [`PB.has_internal_cartesian`](@ref) `true` and implement [`cartesian_to_internal`](@ref), +[`internal_to_cartesian`](@ref), and [`PB.cartesian_size`](@ref). """ PB.AbstractMesh @@ -151,15 +163,8 @@ end Create a CellRange for entire `domain` and supplied `operatorID` """ -function create_default_cellrange(domain::PB.AbstractDomain, grid::Union{PB.AbstractMesh, Nothing}) end +function create_default_cellrange(domain::PB.AbstractDomain, grid::PB.AbstractMeshOrNothing) end -""" - get_region(grid::Union{PB.AbstractMesh, Nothing}, values; selectargs...) -> values_subset, (dim_subset::NamedDimension, ...) - -Return the subset of `values` given by `selectargs` (Grid-specific keywords eg cell=, column=, ...) -and corresponding dimensions (with attached coordinates). -""" -function get_region(grid::Union{PB.AbstractMesh, Nothing}, values) end """ set_subdomain!(grid::PB.AbstractMesh, subdomainname::AbstractString, subdom::PB.AbstractSubdomain, allowcreate::Bool=false) @@ -188,13 +193,9 @@ end # generic handler when subdomainname present -function PB.internal_size(space::Type{<:PB.AbstractSpace}, grid::Union{PB.AbstractMesh, Nothing}, subdomainname::AbstractString) +function PB.internal_size(space::Type{<:PB.AbstractSpace}, grid::PB.AbstractMeshOrNothing, subdomainname::AbstractString) if isempty(subdomainname) - if space === PB.ScalarSpace - return (1, ) - else - return PB.internal_size(space, grid) - end + return PB.internal_size(space, grid) else if space === PB.CellSpace subdomain = get_subdomain(grid, subdomainname) @@ -205,21 +206,134 @@ function PB.internal_size(space::Type{<:PB.AbstractSpace}, grid::Union{PB.Abstra end end +# generic handler when subdomain present +function PB.get_dimensions( + grid::PB.AbstractMeshOrNothing, space::Type{<:PB.AbstractSpace}, subdomainname::AbstractString, + expand_cartesian=false, +) + if isempty(subdomainname) + return PB.get_dimensions(grid, space; expand_cartesian) + else + if space === PB.CellSpace + subdomain = get_subdomain(grid, subdomainname) + return (PB.NamedDimension("subdomain_"*subdomainname, length(subdomain.indices)),) + else + error("get_dimensions: cannot specify subdomain with space=$space (subdomainname=$subdomainname)") + end + end +end + # scalar space always supported -PB.internal_size(space::Type{PB.ScalarSpace}, grid::Union{PB.AbstractMesh, Nothing}) = (1, ) +PB.internal_size(space::Type{PB.ScalarSpace}, grid::PB.AbstractMeshOrNothing) = (1, ) +PB.get_dimensions(grid::PB.AbstractMeshOrNothing, space::Type{PB.ScalarSpace}; expand_cartesian=false) = () # ignore subdomain -PB.internal_size(space::Type{PB.ScalarSpace}, grid::Union{PB.AbstractMesh, Nothing}, subdomainname::AbstractString) = (1, ) +PB.internal_size(space::Type{PB.ScalarSpace}, grid::PB.AbstractMeshOrNothing, subdomainname::AbstractString) = (1, ) + +# fallbacks for cartesian +PB.cartesian_size(space, grid::PB.AbstractMeshOrNothing) = PB.internal_size(space, grid) +PB.has_internal_cartesian(grid::PB.AbstractMeshOrNothing, Space::Type{<:PB.AbstractSpace}) = false + +""" + cartesian_to_internal(grid::PB.AbstractMeshOrNothing, griddata::AbstractArray) + +Map a cartesian (external) array to internal representation +""" +cartesian_to_internal(grid::PB.AbstractMeshOrNothing, griddata::AbstractArray) = griddata + +""" + cartesian_to_internal(grid::PB.AbstractMeshOrNothing, griddata::AbstractArray) + +Map an internal array to cartesian (external) representation +""" +internal_to_cartesian(grid::PB.AbstractMeshOrNothing, griddata::AbstractArray) = griddata + +# generic get all dimensions +function PB.get_dimensions(grid::PB.AbstractMesh) + dims = PB.NamedDimension[] + for space in available_spaces(grid) + append!(dims, PB.get_dimensions(grid, space)) + if PB.has_internal_cartesian(grid, space) + append!(dims, PB.get_dimensions(grid, space; expand_cartesian=true)) + end + end + for sdn in keys(grid.subdomains) + append!(dims, PB.get_dimensions(grid, PB.CellSpace, sdn)) + end + append!(dims, _extra_dimensions(grid)) + return unique(d->d.name, dims) +end + +# implement if subtype has defines additional dimensions not included in list by Space and subdomain +_extra_dimensions(grid::PB.AbstractMesh) = () + +# generic coordinate set / get +function PB.set_coordinates!(grid::PB.AbstractMesh, dimname::AbstractString, coordinates::Vector{<:AbstractString}) + dims = PB.get_dimensions(grid) + idx = findfirst(nd -> nd.name==dimname, dims) + !isnothing(idx) || error("grid does not contain dimension $dimname") + delete!(grid.coordinates, dimname) + if !isempty(coordinates) + grid.coordinates[dimname] = coordinates + end +end + +function PB.get_coordinates(grid::PB.AbstractMesh, dimname::AbstractString) + return get(grid.coordinates, dimname, String[]) +end + +# helper function +function _show_subdomains_dimensions_coordinates(io, grid::PB.AbstractMesh) + println(io, " spaces: ", available_spaces(grid)) + for space in available_spaces(grid) + println(io, " ", space) + if PB.has_internal_cartesian(grid, space) + println(io, " dimensions (internal): ", PB.get_dimensions(grid, space; expand_cartesian=false)) + println(io, " dimensions (cartesian): ", PB.get_dimensions(grid, space; expand_cartesian=true)) + else + println(io, " dimensions: ", PB.get_dimensions(grid, space)) + end + end + println(io, " subdomains: ", keys(grid.subdomains)) + for subdomain in keys(grid.subdomains) + println(io, " dimensions: ", PB.get_dimensions(grid, PB.CellSpace, subdomain)) + end + + extra_dimensions = _extra_dimensions(grid) + if !isempty(extra_dimensions) + println(io, " additional dimensions: ", extra_dimensions) + end + + if !isempty(grid.coordinates) + println(io, " coordinates:") + for (dimname, coordnames) in grid.coordinates + println(io, " ", dimname, " => ", coordnames) + end + end + + return nothing +end + +# generic fallback for named cells +substitute_cell_names(grid::PB.AbstractMesh, cells) = cells +substitute_cell_names(grid::PB.AbstractMesh, cells::Union{Number, Symbol}) = + substitute_cell_names(grid, [cells])[] + +# generic fallback for column indices +function column_indices(grid::PB.AbstractMesh, column) + throw(ArgumentError("grid $grid does not support column selection")) +end + + ################################### # Fallbacks for Domain with no grid ##################################### # allow Vector variables length 1 in a 0D Domain without a grid +available_spaces(grid::Nothing) = (PB.ScalarSpace, PB.CellSpace) PB.internal_size(space::Type{PB.CellSpace}, grid::Nothing) = (1, ) -# allow single dimension -PB.cartesian_size(grid::Nothing) = (1, ) -cartesian_to_internal(grid::Nothing, griddata::AbstractArray) = griddata +PB.get_dimensions(grid::Nothing, space::Type{PB.CellSpace}; expand_cartesian=false) = (PB.NamedDimension("cells", 1), ) get_subdomain(grid::Nothing, subdomainname::AbstractString) = error("get_subdomain: no subdomain $subdomainname") @@ -232,17 +346,7 @@ function create_default_cellrange(domain::PB.AbstractDomain, grid::Nothing; oper return PB.CellRange(domain=domain, indices=1:PB.get_length(domain), operatorID=operatorID) end -""" - get_region(grid::Nothing, values) -> values[] -Fallback for Domain with no grid, assumed 1 cell -""" -function get_region(grid::Nothing, values) - length(values) == 1 || - throw(ArgumentError("grid==Nothing and length(values) != 1")) - return values[], () -end - ################################## # UnstructuredVectorGrid ################################## @@ -258,7 +362,9 @@ Base.@kwdef mutable struct UnstructuredVectorGrid <: PB.AbstractMesh "Define some named cells for plotting (only)" cellnames::Dict{Symbol,Int} = Dict{Symbol,Int}() - subdomains::Dict{String, PB.AbstractSubdomain} = Dict{String, PB.AbstractSubdomain}() + subdomains::Dict{String, PB.AbstractSubdomain} = Dict{String, PB.AbstractSubdomain}() + + coordinates::Dict{String, Vector{String}} = Dict{String, Vector{String}}() end function Base.show(io::IO, grid::UnstructuredVectorGrid) @@ -267,38 +373,19 @@ function Base.show(io::IO, grid::UnstructuredVectorGrid) return nothing end -PB.internal_size(space::Type{PB.CellSpace}, grid::UnstructuredVectorGrid) = (grid.ncells, ) - -# single dimension -PB.cartesian_size(grid::UnstructuredVectorGrid) = (grid.ncells, ) -cartesian_to_internal(grid::UnstructuredVectorGrid, griddata::AbstractArray) = griddata - -""" - get_region(grid::UnstructuredVectorGrid, values; cell=Nothing) -> - values_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `cell::Union{Nothing, Int, Symbol}`: an Int for cell number (first cell = 1), or a Symbol to look up in `cellnames` - `cell = Nothing` is also allowed for a single-cell Domain. -""" -function get_region(grid::UnstructuredVectorGrid, values; cell::Union{Nothing, Int, Symbol}=nothing) - if isnothing(cell) - grid.ncells == 1 || throw(ArgumentError("'cell' argument (an Int or Symbol) is required for an UnstructuredVectorGrid with > 1 cell")) - idx = 1 - elseif cell isa Int - idx = cell - else - idx = get(grid.cellnames, cell, nothing) - !isnothing(idx) || - throw(ArgumentError("cell ':$cell' not present in grid.cellnames=$(grid.cellnames)")) - end - - return ( - values[idx], - (), # no dimensions (ie squeeze out a dimension length 1 for single cell) - ) +function Base.show(io::IO, ::MIME"text/plain", grid::UnstructuredVectorGrid) + println(io, "UnstructuredVectorGrid") + println(io, " ncells=", grid.ncells) + println(io, " cellnames=", grid.cellnames) + _show_subdomains_dimensions_coordinates(io, grid) + return nothing end +available_spaces(grid::UnstructuredVectorGrid) = (PB.ScalarSpace, PB.CellSpace) +PB.internal_size(space::Type{PB.CellSpace}, grid::UnstructuredVectorGrid) = (grid.ncells, ) + +PB.get_dimensions(grid::UnstructuredVectorGrid, space::Type{PB.CellSpace}; expand_cartesian=false) = + (PB.NamedDimension("cells", grid.ncells), ) """ create_default_cellrange(domain, grid::UnstructuredVectorGrid [; operatorID=0]) -> CellRange @@ -309,6 +396,25 @@ function create_default_cellrange(domain::PB.AbstractDomain, grid::UnstructuredV return PB.CellRange(domain=domain, indices=1:grid.ncells, operatorID=operatorID) end +substitute_cell_names(grid::UnstructuredVectorGrid, cells::Union{Number, Symbol}) = + substitute_cell_names(grid, [cells])[] +function substitute_cell_names(grid::UnstructuredVectorGrid, cells) + cell_indices = Int[] + + for cell in cells + if cell isa Int + push!(cell_indices, cell) + else + idx = get(grid.cellnames, cell, nothing) + !isnothing(idx) || + throw(ArgumentError("cell ':$cell' not present in grid.cellnames=$(grid.cellnames)")) + push!(cell_indices, idx) + end + end + + return cell_indices +end + ############################################ # UnstructuredColumnGrid ########################################### @@ -329,12 +435,14 @@ Base.@kwdef mutable struct UnstructuredColumnGrid <: PB.AbstractMesh ncells::Int64 Icolumns::Vector{Vector{Int}} - z_coords::Vector{PB.FixedCoord} = PB.FixedCoord[] + # z_coords::Vector{PB.FixedCoord} = PB.FixedCoord[] "Define optional column names" columnnames::Vector{Symbol} = Symbol[] subdomains::Dict{String, PB.AbstractSubdomain} = Dict{String, PB.AbstractSubdomain}() + + coordinates::Dict{String, Vector{String}} = Dict{String, Vector{String}}() end function Base.show(io::IO, grid::UnstructuredColumnGrid) @@ -343,49 +451,23 @@ function Base.show(io::IO, grid::UnstructuredColumnGrid) return nothing end +function Base.show(io::IO, ::MIME"text/plain", grid::UnstructuredColumnGrid) + println(io, "UnstructuredColumnGrid") + println(io, " ncells=", grid.ncells) + println(io, " columnames=", grid.columnnames) + _show_subdomains_dimensions_coordinates(io, grid) + return nothing +end + +available_spaces(grid::UnstructuredColumnGrid) = (PB.ScalarSpace, PB.ColumnSpace, PB.CellSpace) PB.internal_size(space::Type{PB.CellSpace}, grid::UnstructuredColumnGrid) = (grid.ncells, ) PB.internal_size(space::Type{PB.ColumnSpace}, grid::UnstructuredColumnGrid) = (length(grid.Icolumns), ) -# single dimension -PB.cartesian_size(grid::UnstructuredColumnGrid) = (grid.ncells, ) -cartesian_to_internal(grid::UnstructuredColumnGrid, griddata::AbstractArray) = griddata - -""" - get_region(grid::UnstructuredColumnGrid, values; column, [cell=nothing]) -> - values_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `column::Union{Int, Symbol}`: (may be an Int, or a Symbol to look up in `columnames`) -- `cell::Int`: optional cell index within `column`, highest cell is cell 1 -""" -function get_region(grid::UnstructuredColumnGrid, values; column, cell::Union{Nothing, Int}=nothing) - - if column isa Int - column in 1:length(grid.Icolumns) || - throw(ArgumentError("column index $column out of range")) - colidx = column - else - colidx = findfirst(isequal(column), grid.columnnames) - !isnothing(colidx) || - throw(ArgumentError("columnname '$column' not present in grid.columnnames=$(grid.columnnames)")) - end +PB.get_dimensions(grid::UnstructuredColumnGrid, space::Type{PB.CellSpace}; expand_cartesian=false) = + (PB.NamedDimension("cells", grid.ncells), ) - if isnothing(cell) - indices = grid.Icolumns[colidx] - return ( - values[indices], - (PB.NamedDimension("z", length(indices), PB.get_region(grid.z_coords, indices)), ), - ) - else - # squeeze out z dimension - idx = grid.Icolumns[colidx][cell] - return ( - values[idx], - (), # no dimensions (ie squeeze out a dimension length 1 for single cell) - ) - end - -end +PB.get_dimensions(grid::UnstructuredColumnGrid, space::Type{PB.ColumnSpace}; expand_cartesian=false) = + (PB.NamedDimension("columns", length(grid.Icolumns)), ) """ create_default_cellrange(domain, grid::UnstructuredColumnGrid [; operatorID=0]) -> CellRangeColumns @@ -404,6 +486,20 @@ function create_default_cellrange(domain::PB.AbstractDomain, grid::UnstructuredC ) end +function column_indices(grid::UnstructuredColumnGrid, column) + if column isa Int + column in 1:length(grid.Icolumns) || + throw(ArgumentError("column index $column out of range")) + colidx = column + else + colidx = findfirst(isequal(column), grid.columnnames) + !isnothing(colidx) || + throw(ArgumentError("columnname '$column' not present in grid.columnnames=$(grid.columnnames)")) + end + + return grid.Icolumns[colidx] +end + ###################################### # CartesianLinearGrid ####################################### @@ -430,10 +526,8 @@ Base.@kwdef mutable struct CartesianLinearGrid{N} <: PB.AbstractMesh ncells::Int64 = -1 ncolumns::Int64 = -1 - dimnames::Vector{String} = Vector{String}(undef, N) # netcdf dimension names - dims::Vector{Int} = Vector{Int}(undef, N) - coords::Vector{Vector{Float64}} = Vector{Vector{Float64}}() - coords_edges::Vector{Vector{Float64}} = Vector{Vector{Float64}}() + dimensions::Vector{PB.NamedDimension} = Vector{PB.NamedDimension}(undef, N) # netcdf dimension names + dimensions_extra::Vector{PB.NamedDimension} = [PB.NamedDimension("bnds", 2)] # additional dimension names "index of longitude dimension (if any)" londim::Int = 0 @@ -448,23 +542,47 @@ Base.@kwdef mutable struct CartesianLinearGrid{N} <: PB.AbstractMesh subdomains::Dict{String, PB.AbstractSubdomain} = Dict{String, PB.AbstractSubdomain}() + coordinates::Dict{String, Vector{String}} = Dict{String, Vector{String}}() + "cartesian -> linear index mapping (initialized to `missing` ie linear index contains 0 cells)" - linear_index::Array{Union{Missing,Int32}, N} = Array{Union{Missing,Int32},N}(undef, zeros(Int,N)...) + linear_index::Array{Union{Missing, Int32}, N} = Array{Union{Missing, Int32},N}(undef, zeros(Int,N)...) "linear -> cartesian index mapping (initialized to empty Vector ie linear index contains 0 cells)" cartesian_index::Vector{CartesianIndex{N}} = Vector{CartesianIndex{N}}() end function Base.show(io::IO, grid::CartesianLinearGrid) - dimtuple = NamedTuple{Tuple(Symbol.(grid.dimnames))}(Tuple(grid.dims)) - print(io, "CartesianLinearGrid(ncells=", grid.ncells, ", dimensions: ", dimtuple, + print(io, "CartesianLinearGrid(ncells=", grid.ncells, ", dimensions: ", grid.dimensions, ", dimensions_extra: ", grid.dimensions_extra, ", subdomains: ", keys(grid.subdomains), ")") return nothing end +function Base.show(io::IO, ::MIME"text/plain", grid::CartesianLinearGrid) + println(io, "CartesianLinearGrid") + println(io, " ncells=", grid.ncells) + _show_subdomains_dimensions_coordinates(io, grid) + return nothing +end + +available_spaces(grid::CartesianLinearGrid{2}) = (PB.ScalarSpace, PB.CellSpace) +available_spaces(grid::CartesianLinearGrid{3}) = (PB.ScalarSpace, PB.ColumnSpace, PB.CellSpace) + +PB.has_internal_cartesian(grid::CartesianLinearGrid, ::Type{PB.CellSpace}) = true PB.internal_size(space::Type{PB.CellSpace}, grid::CartesianLinearGrid) = (grid.ncells, ) PB.internal_size(space::Type{PB.ColumnSpace}, grid::CartesianLinearGrid{3}) = (grid.ncolumns, ) +PB.cartesian_size(space::Type{PB.CellSpace}, grid::CartesianLinearGrid) = Tuple(d.size for d in grid.dimensions) + +function PB.get_dimensions(grid::CartesianLinearGrid, space::Type{PB.CellSpace}; expand_cartesian=false) + if expand_cartesian + return grid.dimensions + else + return (PB.NamedDimension("cells", grid.ncells), ) + end +end + +PB.get_dimensions(grid::CartesianLinearGrid{3}, space::Type{PB.ColumnSpace}; expand_cartesian=false) = + (PB.NamedDimension("columns", grid.ncolumns), ) -PB.cartesian_size(grid::CartesianLinearGrid) = Tuple(grid.dims) +_extra_dimensions(grid::CartesianLinearGrid) = grid.dimensions_extra """ @@ -483,10 +601,8 @@ Base.@kwdef mutable struct CartesianArrayGrid{N} <: PB.AbstractMesh ncells::Int64 = -1 - dimnames::Vector{String} = Vector{String}(undef, N) # netcdf dimension names - dims::Vector{Int} = Vector{Int}(undef, N) - coords::Vector{Vector{Float64}} = Vector{Vector{Float64}}() - coords_edges::Vector{Vector{Float64}} = Vector{Vector{Float64}}() + dimensions::Vector{PB.NamedDimension} = Vector{PB.NamedDimension}(undef, N) # netcdf dimension names + dimensions_extra::Vector{PB.NamedDimension} = [PB.NamedDimension("bnds", 2)] # additional dimension names "index of longitude dimension (if any)" londim::Int = 0 @@ -500,108 +616,50 @@ Base.@kwdef mutable struct CartesianArrayGrid{N} <: PB.AbstractMesh display_mult::Vector{Float64} = ones(N) subdomains::Dict{String, PB.AbstractSubdomain} = Dict{String, PB.AbstractSubdomain}() + + coordinates::Dict{String, Vector{String}} = Dict{String, Vector{String}}() end function Base.show(io::IO, grid::CartesianArrayGrid) - dimtuple = NamedTuple{Tuple(Symbol.(grid.dimnames))}(Tuple(grid.dims)) - print(io, "CartesianArrayGrid(ncells=", grid.ncells, ", dimensions: ", dimtuple, + print(io, "CartesianArrayGrid(ncells=", grid.ncells, ", dimensions: ", grid.dimensions, ", dimensions_extra ", grid.dimensions_extra, ", subdomains: ", keys(grid.subdomains), ")") return nothing end -PB.internal_size(space::Type{PB.CellSpace}, grid::CartesianArrayGrid) = Tuple(grid.dims) - -PB.cartesian_size(grid::CartesianArrayGrid) = Tuple(grid.dims) - - -""" - get_region(grid::Union{CartesianLinearGrid{2}, CartesianArrayGrid{2}} , internalvalues; [i=i_idx], [j=j_idx]) -> - arrayvalues_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `i::Int`: optional, slice along first dimension -- `j::Int`: optional, slice along second dimension - -`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to -an Array (2D if neither i, j arguments present, 1D if i or j present, 0D ie one cell if both present) -""" -function get_region( - grid::Union{CartesianLinearGrid{2}, CartesianArrayGrid{2}}, internalvalues; - i::Union{Integer, Colon}=Colon(), j::Union{Integer, Colon}=Colon() -) - return _get_region(grid, internalvalues, [i, j]) -end - -""" - get_region(grid::Union{CartesianLinearGrid{3}, CartesianArrayGrid{3}}, internalvalues; [i=i_idx], [j=j_idx]) -> - arrayvalues_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `i::Int`: optional, slice along first dimension -- `j::Int`: optional, slice along second dimension -- `k::Int`: optional, slice along third dimension - -`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to -an Array (3D if neither i, j, k arguments present, 2D if one of i, j or k present, 1D if two present, -0D ie one cell if i, j, k all specified). -""" -function get_region( - grid::Union{CartesianLinearGrid{3}, CartesianArrayGrid{3}}, internalvalues; - i::Union{Integer, Colon}=Colon(), j::Union{Integer, Colon}=Colon(), k::Union{Integer, Colon}=Colon() -) - return _get_region(grid, internalvalues, [i, j, k]) +function Base.show(io::IO, ::MIME"text/plain", grid::CartesianArrayGrid) + println(io, "CartesianArrayGrid") + println(io, " ncells=", grid.ncells) + _show_subdomains_dimensions_coordinates(io, grid) + return nothing end -function _get_region( - grid::Union{CartesianLinearGrid, CartesianArrayGrid}, internalvalues, indices -) - if !isempty(grid.coords) && !isempty(grid.coords_edges) - dims = [ - PB.NamedDimension(grid.dimnames[idx], grid.coords[idx], grid.coords_edges[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - elseif !isempty(grid.coords) - dims = [ - PB.NamedDimension(grid.dimnames[idx], grid.coords[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - else - dims = [ - PB.NamedDimension(grid.dimnames[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - end +available_spaces(grid::CartesianArrayGrid) = (PB.ScalarSpace, PB.CellSpace) +PB.internal_size(space::Type{PB.CellSpace}, grid::CartesianArrayGrid) = Tuple(d.size for d in grid.dimensions) - values = internal_to_cartesian(grid, internalvalues) - if !all(isequal(Colon()), indices) - values = values[indices...] - end +PB.get_dimensions(grid::CartesianArrayGrid, space::Type{PB.CellSpace}; expand_cartesian=false) = + grid.dimensions - return values, Tuple(dims) -end +_extra_dimensions(grid::CartesianArrayGrid) = grid.dimensions_extra """ CartesianGrid( - gridtype, dimnames::Vector{<:AbstractString}, dims, coords, coords_edges; + gridtype, dimensions::Vector{PB.NamedDimension}; [londim] [, latdim] [,zdim=0], [,zidxsurface=0], [, ztoheight=1.0]) -> grid::CartesianLinearGrid Create a CartesianLinearGrid or CartesianArrayGrid from dimensions and coordinates. """ function CartesianGrid( GridType::Type{<:Union{CartesianLinearGrid, CartesianArrayGrid}}, - dimnames::Vector{<:AbstractString}, dims, coords=Vector{Vector{Float64}}(), coords_edges=Vector{Vector{Float64}}(); - londim=findfirst(isequal("lon"), dimnames), - latdim=findfirst(isequal("lat"), dimnames), + dimensions::Vector{PB.NamedDimension}; + londim=findfirst(d-> d.name == "lon", dimensions), + latdim=findfirst(d-> d.name == "lat", dimensions), zdim=0, zidxsurface=0, ztoheight=1.0 ) - grid = GridType{length(dimnames)}( - dimnames=dimnames, - dims=dims, - coords=coords, - coords_edges=coords_edges, + grid = GridType{length(dimensions)}(; + dimensions=dimensions, londim=londim, latdim=latdim, zdim=zdim, @@ -612,21 +670,24 @@ function CartesianGrid( grid.display_mult[zdim] = ztoheight end + sz = Tuple(d.size for d in dimensions) if GridType == CartesianLinearGrid # create linear index (with no points) - grid.linear_index = Array{Union{Missing,Int32},length(grid.dims)}(undef, grid.dims...) + grid.linear_index = Array{Union{Missing, Int32},length(grid.dimensions)}(undef, sz...) grid.ncells = 0 else - grid.ncells = prod(grid.dims) + grid.ncells = prod(sz) end return grid end """ - CartesianGrid(griddtype, ncfilename::AbstractString, dimnames::Vector{<:AbstractString}; equalspacededges=false) -> grid::CartesianLinearGrid + CartesianGrid(GridType, ncfilename::AbstractString, dimnames::Vector{<:AbstractString}; equalspacededges=false) - + > grid::GridType, coord_vars, coord_edges_vars Read a netcdf file with CF1.0 coordinates, and create corresponding a CartesianLinearGrid or CartesianArrayGrid from dimensions `dimnames`. +`coord_vars` and optional `coord_edges_vars` are Vectors of `coorddimname=>values` for creating coordinate variables. """ function CartesianGrid( GridType::Type{<:Union{CartesianLinearGrid, CartesianArrayGrid}}, @@ -637,29 +698,36 @@ function CartesianGrid( @info "CartesianGrid creating $GridType{$(length(dimnames))} from dimnames=$dimnames in netcdf file $(ncfilename)" grid = GridType{length(dimnames)}() - grid.coords = Vector{Vector{Float64}}(undef, length(dimnames)) - grid.coords_edges = Vector{Vector{Float64}}(undef, length(dimnames)) + + coord_vars = [] + coord_edges_vars = [] NCDatasets.Dataset(ncfilename) do ds # read dimensions and coordinates for i in eachindex(dimnames) - dimname = dimnames[i] - grid.dimnames[i] = dimname - grid.dims[i] = ds.dim[dimname] - @info " read dimension $(dimname) = $(grid.dims[i])" + dimname = dimnames[i] + dim = ds.dim[dimname] + @info " read dimension $dimname = $dim" + grid.dimensions[i] = PB.NamedDimension(dimname, dim) v_cf = ds[dimname] - grid.coords[i] = Array(v_cf) + coord_vals = Array(v_cf) + push!(coord_vars, dimname=>coord_vals) if haskey(v_cf.attrib, "edges") edgesname = v_cf.attrib["edges"] - @info " reading coordinate edges from $edgesname" - grid.coords_edges[i] = Array(ds[edgesname]) + edgesdim = ds.dim[edgesname] + @info " reading coordinate edges from $edgesname = $edgesdim" + push!(grid.dimensions_extra, PB.NamedDimension(edgesname, edgesdim)) + push!(coord_edges_vars, edgesname=>Array(ds[edgesname])) elseif equalspacededges - es = grid.coords[i][2] - grid.coords[i][1] - @info " assuming equal spacing $es to calculate coordinate edges" - grid.coords_edges[i] = [grid.coords[i] .- es/2.0; grid.coords[i][end] + es/2.0 ] + edgesname = dimname*"_edges" + ew = coord_vals[2] - coord_vals[1] + edgesvals = [coord_vals .- ew/2.0; coord_vals[end] + ew/2.0 ] + @info " assuming equal spacing $ew to calculate coordinate edges $edgesname" + push!(grid.dimensions_extra, PB.NamedDimension(edgesname, dim+1)) + push!(coord_edges_vars, edgesname=>edgesvals) else - error(" no edges attribute and equalspacededges=false") + @info " no edges attribute and equalspacededges=false" end if v_cf.attrib["standard_name"] == "longitude" @@ -672,25 +740,26 @@ function CartesianGrid( @info " dim $i got standard_name=='depth', using this as z dimension" grid.zdim = i grid.display_mult[i] = -1.0 - if grid.coords[i][1] < grid.coords[i][end] + if coord_vals[1] < coord_vals[end] grid.zidxsurface = 1 else - grid.zidxsurface = length(grid.coords[i]) + grid.zidxsurface = length(coord_vals) end @info " surface is index $(grid.zidxsurface)" end end end + grid_size = Tuple(nd.size for nd in grid.dimensions) if GridType == CartesianLinearGrid - # create linear index (with no points) - grid.linear_index = Array{Union{Missing,Int32},length(grid.dims)}(undef, grid.dims...) + # create linear index (with no points) + grid.linear_index = Array{Union{Missing,Int32}, length(grid_size)}(undef, grid_size...) grid.ncells = 0 else - grid.ncells = prod(grid.dims) + grid.ncells = prod(grid_size) end - return grid + return grid, coord_vars, coord_edges_vars end """ @@ -703,7 +772,8 @@ function set_linear_index(grid::CartesianLinearGrid{3}, v_i, v_j, v_k) grid.ncells = length(v_i) grid.ncolumns = 0 - grid.linear_index = Array{Union{Missing,Int32}, 3}(undef, grid.dims...) + grid_size = Tuple(nd.size for nd in grid.dimensions) + grid.linear_index = Array{Union{Missing, Int32}, 3}(undef, grid_size...) grid.cartesian_index = Vector{CartesianIndex{3}}(undef, grid.ncells) fill!(grid.linear_index, missing) for l in eachindex(v_i) @@ -726,7 +796,8 @@ function set_linear_index(grid::CartesianLinearGrid{2}, v_i, v_j) grid.ncells = length(v_i) - grid.linear_index = Array{Union{Missing,Int32}, 2}(undef, grid.dims...) + grid_size = Tuple(nd.size for nd in grid.dimensions) + grid.linear_index = Array{Union{Missing,Int32}, 2}(undef, grid_size) grid.cartesian_index = Vector{CartesianIndex{2}}(undef, grid.ncells) fill!(grid.linear_index, missing) for l in eachindex(v_i) @@ -750,76 +821,55 @@ function cartesian_to_internal(grid::CartesianLinearGrid, griddata::AbstractArra return [griddata[i] for i in grid.cartesian_index] end -cartesian_to_internal(grid::CartesianArrayGrid, griddata::AbstractArray) = griddata - """ - internal_to_cartesian(grid::CartesianLinearGrid, internaldata::AbstractVector [,missing_value=missing]) -> griddata::Array + internal_to_cartesian(grid::CartesianLinearGrid, internaldata::AbstractArray [,missing_value=missing]) -> griddata::Array - Convert Vector `internaldata` (on `grid.linear_index`) to Cartesian Array `griddata` (with `missing_value` where no data) + Convert Array `internaldata` (first index on `grid.linear_index`, up to 2 additional indices for additional data dimensions) + to Cartesian Array `griddata` (with `missing_value` where no data) """ -function internal_to_cartesian(grid::CartesianLinearGrid, internaldata::AbstractArray; missing_value=missing) - grid.ncells == length(internaldata) || error("grid and data size mismatch") +function internal_to_cartesian(grid::CartesianLinearGrid, internaldata::AbstractArray{T, N}; missing_value=missing) where {T, N} + grid.ncells == size(internaldata, 1) || error("grid and data size mismatch") if missing_value isa Missing - gr_eltype = Union{Missing, eltype(internaldata)} + gr_eltype = Union{Missing, T} else - gr_eltype = eltype(internaldata) + gr_eltype = T + end + + if N == 1 # internaldata is a Vector; no extra data dimensions + griddata = similar(grid.linear_index, gr_eltype) + fill!(griddata, missing_value) + griddata[grid.cartesian_index] .= internaldata + elseif N == 2 # 1 extra data dimension`` + griddata = similar(grid.linear_index, gr_eltype, size(grid.linear_index)..., size(internaldata, 2)) + fill!(griddata, missing_value) + griddata[grid.cartesian_index, :] .= internaldata + elseif N == 3 # 2 extra data dimensions + griddata = similar(grid.linear_index, gr_eltype, size(grid.linear_index)..., size(internaldata, 2), size(internaldata, 3)) + fill!(griddata, missing_value) + griddata[grid.cartesian_index, :, :] .= internaldata + else + error("data with $(N-1) extra data dimensions not supported") end - griddata = similar(grid.linear_index, gr_eltype) - fill!(griddata, missing_value) - griddata[grid.cartesian_index] .= internaldata return griddata end -internal_to_cartesian(grid::CartesianArrayGrid, internaldata::AbstractArray; missing_value=missing) = internaldata -function get_lon(grid::CartesianLinearGrid, linear_idx::Integer) +function get_lon_idx(grid::CartesianLinearGrid, linear_idx::Integer) cartesian_idx = grid.cartesian_index[linear_idx] - return _get_coord(grid.coords, grid.londim, "longitude", cartesian_idx) + return cartesian_idx[grid.londim] end -get_lon(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = - _get_coord(grid.coords, grid.londim, "longitude", cartesian_idx) +get_lon_idx(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = cartesian_idx[grid.londim] -function get_lon_edges(grid::CartesianLinearGrid, linear_idx::Integer) +function get_lat_idx(grid::CartesianLinearGrid, linear_idx::Integer) cartesian_idx = grid.cartesian_index[linear_idx] - return _get_coord_edges(grid.coords_edges, grid.londim, "longitude", cartesian_idx) + return cartesian_idx[grid.latdim] end -get_lon_edges(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = - _get_coord_edges(grid.coords_edges, grid.londim, "longitude", cartesian_idx) - - -function get_lat(grid::CartesianLinearGrid, linear_idx::Integer) - cartesian_idx = grid.cartesian_index[linear_idx] - return _get_coord(grid.coords, grid.latdim, "latitude", cartesian_idx) -end +get_lat_idx(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = cartesian_idx[grid.latdim] -get_lat(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = - _get_coord(grid.coords, grid.latdim, "latitude", cartesian_idx) - -function get_lat_edges(grid::CartesianLinearGrid, linear_idx::Integer) - cartesian_idx = grid.cartesian_index[linear_idx] - return _get_coord_edges(grid.coords_edges, grid.latdim, "latitude", cartesian_idx) -end - -get_lat_edges(grid::CartesianArrayGrid{N}, cartesian_idx::CartesianIndex{N}) where {N} = - _get_coord_edges(grid.coords_edges, grid.latdim, "latitude", cartesian_idx) - - -function _get_coord(coords, cdim, dimname, cartesian_idx::CartesianIndex) - cdim > 0 || error("grid has no $dimname dimension") - return coords[cdim][cartesian_idx[cdim]] -end - -function _get_coord_edges(coords_edges, cdim, dimname, cartesian_idx::CartesianIndex) - cdim > 0 || error("grid has no $dimname dimension") - return ( - coords_edges[cdim][cartesian_idx[cdim]], - coords_edges[cdim][cartesian_idx[cdim]+1], - ) -end """ linear_indices_cartesian_ranges(grid::CartesianLinearGrid, rangestuple) -> lindices @@ -851,7 +901,7 @@ Create a range of cells within a region of CartesianLinearGrid specified by `ran """ function cellrange_cartesiantile(domain, grid::CartesianLinearGrid{2}, rangestuple; operatorID=0) indices = Int[] - r1, r2 = _expandrangetuple(grid.dims, rangestuple[1:2]) + r1, r2 = _expandrangetuple(grid.dimensions, rangestuple[1:2]) for i=1:length(grid.cartesian_index) ci = grid.cartesian_index[i] @@ -870,9 +920,9 @@ function cellrange_cartesiantile(domain, grid::CartesianLinearGrid{3}, rangestup # recreate (all) surface indices surfindices = [ci for ci in grid.cartesian_index if ci[3] == grid.zidxsurface] - r1, r2, r3 = _expandrangetuple(grid.dims, rangestuple) + r1, r2, r3 = _expandrangetuple(grid.dimensions, rangestuple) indices = Int[] - colindices = Vector{Pair{Int,Vector{Int}}}() + colindices = Vector{Pair{Int, Vector{Int}}}() # iterate through surface indices and generate columns for is in 1:length(surfindices) sci = surfindices[is] @@ -908,13 +958,13 @@ function cellrange_cartesiantile(domain, grid::CartesianLinearGrid{3}, rangestup end function cellrange_cartesiantile(domain, grid::CartesianArrayGrid{2}, rangestuple; operatorID=0) - rt = _expandrangetuple(grid.dims, rangestuple[1:2]) + rt = _expandrangetuple(grid.dimensions, rangestuple[1:2]) return PB.CellRange(domain=domain, indices=CartesianIndices(rt), operatorID=operatorID) end function cellrange_cartesiantile(domain, grid::CartesianArrayGrid{3}, rangestuple; operatorID=0) - rt = _expandrangetuple(grid.dims, rangestuple[1:3]) + rt = _expandrangetuple(grid.dimensions, rangestuple[1:3]) return PB.CellRange(domain=domain, indices=CartesianIndices(rt), operatorID=operatorID) end @@ -929,12 +979,12 @@ function create_default_cellrange(domain::PB.AbstractDomain, grid::Union{Cartesi end "expand (1:2, :, :) replacing : with ranges" -function _expandrangetuple(dims, rangestuple) +function _expandrangetuple(nameddims::Vector{PB.NamedDimension}, rangestuple) ranges = [] for i = 1:length(rangestuple) r = rangestuple[i] if r isa Colon - push!(ranges, 1:dims[i]) + push!(ranges, 1:nameddims[i].size) else push!(ranges, r) end @@ -1006,10 +1056,11 @@ traversal of the 3D grid in order k, j, i """ function linear_kji_order(grid, v_i, v_j, v_k) - length(grid.dims) == 3 || error("linear_kji_order grid is not a 3D grid") + length(grid.dimensions) == 3 || error("linear_kji_order grid is not a 3D grid") # construct a 3D Array and fill with linear index values - linear_index = Array{Union{Missing,Int32},length(grid.dims)}(missing, grid.dims...) + grid_size = Tuple(nd.size for nd in grid.dimensions) + linear_index = Array{Union{Missing,Int32},length(grid_size)}(missing, grid_size...) for l in eachindex(v_i) linear_index[v_i[l], v_j[l], v_k[l]] = l end @@ -1017,9 +1068,9 @@ function linear_kji_order(grid, v_i, v_j, v_k) # visit 3D Array in k, j, i order and record linear index values perm = Vector{Int}(undef, length(v_i)) ip = 0 - for i in 1:grid.dims[1] - for j in 1:grid.dims[2] - for k in 1:grid.dims[3] + for i in 1:grid.dimensions[1].size + for j in 1:grid.dimensions[2].size + for k in 1:grid.dimensions[3].size if !ismissing(linear_index[i, j, k]) ip += 1 perm[ip] = linear_index[i, j, k] diff --git a/src/Model.jl b/src/Model.jl index 29b35e07..4c0139dc 100644 --- a/src/Model.jl +++ b/src/Model.jl @@ -924,9 +924,18 @@ end # multiline form function Base.show(io::IO, ::MIME"text/plain", model::Model) println(io, "Model") - println(io, "\tname='", model.name,"'") - println(io, "\tconfig_files='", model.config_files,"'") - println(io, "\tdomains=", model.domains) + println(io, " name: '", model.name,"'") + println(io, " config_files: '", model.config_files,"'") + println(io, " domains:") + for dom in model.domains + iodom = IOBuffer() + show(iodom, MIME"text/plain"(), dom; show_reactions=false, show_variables=false) + seekstart(iodom) + for line in eachline(iodom) + println(io, " ", line) + end + end + return nothing end """ diff --git a/src/PALEOboxes.jl b/src/PALEOboxes.jl index 8bad61c9..09992e3a 100644 --- a/src/PALEOboxes.jl +++ b/src/PALEOboxes.jl @@ -203,7 +203,7 @@ end rdict = find_all_reactions() reactionlist = [ "ReactionFluxTransfer", "ReactionReservoirScalar", "ReactionFluxPerturb", "ReactionReservoir", - "ReactionReservoirForced", "ReactionSum", "ReactionFluxTarget", "ReactionForceInterp", "ReactionGrid2DNetCDF", + "ReactionReservoirForced", "ReactionSum", "ReactionFluxTarget", "ReactionForceInterp", "ReactionAreaVolumeValInRange", "ReactionReservoirWellMixed", "ReactionForceGrid", "ReactionConst", "ReactionRestore", "ReactionScalarConst", "ReactionVectorSum", "ReactionWeightedMean", "ReactionReservoirTotal", "ReactionUnstructuredVectorGrid", "ReactionCartesianGrid", "ReactionReservoirConst", diff --git a/src/Types.jl b/src/Types.jl index eae29f5e..94cc1c22 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -1,6 +1,17 @@ # TODO there doesn't seem to be an easy way of parameterising ModelData by an Array type ? const PaleoArrayType = Array +################################ +# CommonDataModel adaptor +################################ + +""" + CDModel(x) + +Create a CommonDataModel adaptor for PALEO object x +""" +function CDModel end + ################################ # Parameters ############################### @@ -100,6 +111,8 @@ end abstract type AbstractMesh end +const AbstractMeshOrNothing = Union{AbstractMesh, Nothing} + """ function get_mesh(obj, ...) @@ -155,6 +168,7 @@ Return a Tables.jl data table view for PALEO object `obj` """ function get_table end + ############################################# # ReactionMethodDispatchList ############################################# diff --git a/src/VariableDomain.jl b/src/VariableDomain.jl index b20f1b6c..287c06be 100644 --- a/src/VariableDomain.jl +++ b/src/VariableDomain.jl @@ -101,6 +101,14 @@ function host_dependent(var::VariableDomContribTarget) return isnothing(var.var_target) end +function get_dimensions(v::VariableDomain; expand_cartesian=false) + dims = NamedDimension[] + append!(dims, PB.get_dimensions(v.domain.grid, v.attributes[:space])) + for ddn in v.attributes[:data_dims] + append!(get_data_dimension(v.domain, ddn)) + end + return dims +end ################################################ # Field and data access @@ -123,15 +131,13 @@ end """ set_data!(var::VariableDomain, modeldata, arrays_idx::Int, data) -Set [`VariableDomain`](@ref) to a Field containing `data`. - -Calls [`wrap_field`](@ref) to create a new [`Field`](@ref). +Set [`VariableDomain`](@ref) to a newly-created [`Field`](@ref) containing `data`. """ function set_data!(var::VariableDomain, modeldata::AbstractModelData, arrays_idx::Int, data) variable_data = get_domaindata(modeldata, var.domain, arrays_idx).variable_data - variable_data[var.ID] = wrap_field( + variable_data[var.ID] = Field( data, get_attribute(var, :field_data), Tuple(get_data_dimension(domain, dn) for dn in get_attribute(var, :data_dims)), @@ -253,7 +259,7 @@ function allocate_variables!( check_units(v; check_units_opt) data_dims = Tuple( - get_data_dimension(v.domain, dimname) + get_dimension(v.domain, dimname) for dimname in get_attribute(v, :data_dims) ) @@ -290,7 +296,7 @@ function allocate_variables!( v_field = get_field(v, modeldata, 1) else # allocate new field array - v_field = allocate_field( + v_field = Field( field_data, data_dims, mdeltype, space, v.domain.grid; allocatenans=modeldata.allocatenans ) diff --git a/src/reactioncatalog/GridForcings.jl b/src/reactioncatalog/GridForcings.jl index 7e798556..90d0be56 100644 --- a/src/reactioncatalog/GridForcings.jl +++ b/src/reactioncatalog/GridForcings.jl @@ -170,7 +170,7 @@ function _prepare_data(rj::ReactionForceGrid, ds) # NB: A PALEO Cartesian Grid may define an internal_size for model Variables that has different dimensions (eg a linear Vector) # to the n-D cartesian_size of the forcings read from the NetCDF file - ncartesiandims = length(PB.cartesian_size(rj.domain.grid)) # as read from NetCDF + ncartesiandims = length(PB.cartesian_size(PB.CellSpace, rj.domain.grid)) # as read from NetCDF cartesiancolons = fill(Colon(), ncartesiandims) ninternaldims = length(PB.internal_size(PB.CellSpace, rj.domain.grid)) # PALEO array layout internalcolons = fill(Colon(), ninternaldims) diff --git a/src/reactioncatalog/GridReactions.jl b/src/reactioncatalog/GridReactions.jl index c022c0fa..6a73d383 100644 --- a/src/reactioncatalog/GridReactions.jl +++ b/src/reactioncatalog/GridReactions.jl @@ -63,7 +63,13 @@ function PB.set_model_geometry(rj::ReactionCartesianGrid, model::PB.Model) grid = PB.Grids.CartesianGrid( PB.Grids.CartesianArrayGrid, - rj.pars.dimnames.v, rj.pars.dims.v + [ + PB.NamedDimension(name, sz) + for (name, sz) in PB.IteratorUtils.zipstrict( + rj.pars.dimnames.v, rj.pars.dims.v; + errmsg="parameters 'dimnames' and 'dims' are not the same length" + ) + ] ) rj.domain.grid = grid @@ -111,13 +117,15 @@ Base.@kwdef mutable struct ReactionGrid2DNetCDF{P} <: PB.AbstractReaction description="radius to calculate cell area from spherical geometry (if area_var = \"\")"), ) + coord_varnames_values = nothing + coord_edges_varnames_values = nothing end function PB.set_model_geometry(rj::ReactionGrid2DNetCDF, model::PB.Model) @info "set_model_geometry $(PB.fullname(rj))" @info " reading 2D grid from $(rj.pars.grid_file[])" - grid2D = PB.Grids.CartesianGrid( + grid2D, rj.coord_varnames_values, rj.coord_edges_varnames_values = PB.Grids.CartesianGrid( rj.pars.grid_type[], rj.pars.grid_file[], rj.pars.coordinate_names.v, equalspacededges=rj.pars.equalspacededges[] @@ -125,8 +133,8 @@ function PB.set_model_geometry(rj::ReactionGrid2DNetCDF, model::PB.Model) if rj.pars.grid_type[] == PB.Grids.CartesianLinearGrid # define a linear index including every cell, in column-major order (first indices are consecutive in memory) - v_i = vec([i for i=1:grid2D.dims[1], j=1:grid2D.dims[2]]) - v_j = vec([j for i=1:grid2D.dims[1], j=1:grid2D.dims[2]]) + v_i = vec([i for i=1:grid2D.dimensions[1].size, j=1:grid2D.dimensions[2].size]) + v_j = vec([j for i=1:grid2D.dimensions[1].size, j=1:grid2D.dimensions[2].size]) PB.Grids.set_linear_index(grid2D, v_i, v_j) end @@ -145,10 +153,22 @@ function PB.register_methods!(rj::ReactionGrid2DNetCDF) PB.VarPropScalarStateIndep("Asurf_total", "m^2", "total horizontal area of surface"), ] + coord_vars = [ + PB.VarPropScalarStateIndep(coord_name, "", "coordinate variable"; + attributes=(:field_data=>PB.ArrayScalarData, :data_dims=>(coord_name,), )) + for (coord_name, _) in rj.coord_varnames_values + ] + + coord_edges_vars = [ + PB.VarPropScalarStateIndep(coord_name, "", "coordinate edge variable"; + attributes=(:field_data=>PB.ArrayScalarData, :data_dims=>(coord_name,), )) + for (coord_name, _) in rj.coord_edges_varnames_values + ] + PB.add_method_setup!( rj, setup_grid_2DNetCDF, - (PB.VarList_namedtuple(grid_vars),) + (PB.VarList_namedtuple(grid_vars), PB.VarList_vector(coord_vars), PB.VarList_vector(coord_edges_vars)) ) return nothing @@ -157,7 +177,7 @@ end function setup_grid_2DNetCDF( m::PB.ReactionMethod, pars, - (grid_vars, ), + (grid_vars, coord_vars, coord_edges_vars), cellrange::PB.AbstractCellRange, attribute_name ) @@ -172,20 +192,32 @@ function setup_grid_2DNetCDF( length(cellrange.indices) == grid2D.ncells || error("tiled cellrange not supported") + for (coord_var, (coord_name, coord_values)) in PB.IteratorUtils.zipstrict(coord_vars, rj.coord_varnames_values) + coord_var .= coord_values + end + + for (coord_edge_var, (coord_edge_name, coord_edge_values)) in PB.IteratorUtils.zipstrict(coord_edges_vars, rj.coord_edges_varnames_values) + coord_edge_var .= coord_edge_values + end + if !isempty(pars.area_var[]) NCDatasets.Dataset(pars.grid_file[]) do ds area2D = Array(ds[pars.area_var[]][:, :, 1]) - grid_vars.Asurf .= cartesian_to_internal(rj.domain.grid, area2D) + grid_vars.Asurf .= PB.Grids.cartesian_to_internal(rj.domain.grid, area2D) end - else - for i in cellrange.indices - lon_edges = PB.Grids.get_lon_edges(grid2D, i) - lat_edges = PB.Grids.get_lat_edges(grid2D, i) + elseif !isempty(coord_edges_vars) + lonedgevar, latedgevar = coord_edges_vars[[grid2D.londim, grid2D.latdim]] + for idx in cellrange.indices + lonidx = PB.Grids.get_lon_idx(grid2D, idx) + latidx = PB.Grids.get_lat_idx(grid2D, idx) + lon_edges = lonedgevar[lonidx:lonidx+1] + lat_edges = latedgevar[latidx:latidx+1] area = calc_spherical_area(pars.planet_radius[], lon_edges, lat_edges) - grid_vars.Asurf[i] = area + grid_vars.Asurf[idx] = area end + else + @warn "no area_var or coords edges specified, not calculating Asurf" end - grid_vars.Asurf_total[] = sum(grid_vars.Asurf) diff --git a/test/runfieldtests.jl b/test/runfieldtests.jl index e4f44cdc..10323a0f 100644 --- a/test/runfieldtests.jl +++ b/test/runfieldtests.jl @@ -9,7 +9,7 @@ import PALEOboxes as PB @testset "ScalarData, ScalarSpace" begin - f = PB.allocate_field(PB.ScalarData, (), Float64, PB.ScalarSpace, nothing; allocatenans=true) + f = PB.Field(PB.ScalarData, (), Float64, PB.ScalarSpace, nothing; allocatenans=true) @test isnan(f.values[]) @@ -18,7 +18,7 @@ import PALEOboxes as PB @testset "ScalarData, CellSpace, no grid" begin # check that a CellSpace Field with no grid behaves as a ScalarSpace Field - f = PB.allocate_field(PB.ScalarData, (), Float64, PB.CellSpace, nothing; allocatenans=true) + f = PB.Field(PB.ScalarData, (), Float64, PB.CellSpace, nothing; allocatenans=true) @test isnan(f.values[]) end @@ -26,14 +26,14 @@ import PALEOboxes as PB @testset "ScalarData, CellSpace, UnstructuredColumnGrid" begin g = PB.Grids.UnstructuredColumnGrid(ncells=5, Icolumns=[collect(1:5)]) - f = PB.allocate_field(PB.ScalarData, (), Float64, PB.CellSpace, g; allocatenans=true) + f = PB.Field(PB.ScalarData, (), Float64, PB.CellSpace, g; allocatenans=true) @test isnan.(f.values) == [true for i in 1:5] end @testset "ArrayScalarData, ScalarSpace" begin - f = PB.allocate_field(PB.ArrayScalarData, (PB.NamedDimension("test", 2, []), ), Float64, PB.ScalarSpace, nothing; allocatenans=true) + f = PB.Field(PB.ArrayScalarData, (PB.NamedDimension("test", 2), ), Float64, PB.ScalarSpace, nothing; allocatenans=true) @test isnan.(f.values) == [true, true] @@ -42,7 +42,7 @@ import PALEOboxes as PB @testset "ArrayScalarData, CellSpace, no grid" begin # TODO this is possibly inconsistent with (ScalarData, CellSpace, no grid), # as Field.values here is a (1, 2) Array, not a (2,) Vector - f = PB.allocate_field(PB.ArrayScalarData, (PB.NamedDimension("test", 2, []), ), Float64, PB.CellSpace, nothing; allocatenans=true) + f = PB.Field(PB.ArrayScalarData, (PB.NamedDimension("test", 2), ), Float64, PB.CellSpace, nothing; allocatenans=true) @test_broken size(f.values) == (2, ) # TODO should be a Vector ? @test size(f.values) == (1, 2) @@ -53,7 +53,7 @@ import PALEOboxes as PB @testset "ArrayScalarData, CellSpace, UnstructuredColumnGrid" begin g = PB.Grids.UnstructuredColumnGrid(ncells=5, Icolumns=[collect(1:5)]) - f = PB.allocate_field(PB.ArrayScalarData, (PB.NamedDimension("test", 2, []), ), Float64, PB.CellSpace, g; allocatenans=true) + f = PB.Field(PB.ArrayScalarData, (PB.NamedDimension("test", 2), ), Float64, PB.CellSpace, g; allocatenans=true) @test size(f.values) == (5, 2) From b97581c8aa7fe53a46284533929bf841719ac31a Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Sun, 29 Dec 2024 22:41:59 +0000 Subject: [PATCH 2/6] Tidy up, fix doc build --- docs/src/DomainsVariablesFields.md | 40 +++++++++++++++++------------- src/CoordsDims.jl | 13 +++++----- src/Grids.jl | 29 +++++++++++----------- 3 files changed, 45 insertions(+), 37 deletions(-) diff --git a/docs/src/DomainsVariablesFields.md b/docs/src/DomainsVariablesFields.md index 72503bba..aa811d58 100644 --- a/docs/src/DomainsVariablesFields.md +++ b/docs/src/DomainsVariablesFields.md @@ -21,6 +21,25 @@ CurrentModule = PALEOboxes Domain ``` +### Dimensions and Coordinates + +PALEO approximately follows the [Common Data Model](https://docs.unidata.ucar.edu/netcdf-c/current/netcdf_data_model.html) used by NetCDF and other geoscience data formats. + +Domains provide named dimensions, to which coordinate variables may be attached (these are just normal PALEO variables). +- Domain spatial dimensions are provided by the Domain grid (see `Grids`) +- Additional Domain data dimensions (eg a wavelength grid) may be set by Reactions (see [`set_data_dimension!`](@ref)). + +```@meta +CurrentModule = PALEOboxes +``` +```@docs +NamedDimension +get_dimensions +get_dimension +get_coordinates +set_coordinates! +``` + ### Grids ```@meta CurrentModule = PALEOboxes @@ -28,8 +47,12 @@ CurrentModule = PALEOboxes ```@docs AbstractMesh +Grids.available_spaces +has_internal_cartesian internal_size cartesian_size +Grids.cartesian_to_internal +Grids.internal_to_cartesian ``` ```@meta CurrentModule = PALEOboxes.Grids @@ -61,23 +84,6 @@ subdomain_view subdomain_indices ``` -#### Regions, Dimensions and Coordinates -```@meta -CurrentModule = PALEOboxes.Grids -``` -```@docs -get_region -``` - -Grids may define name dimensions and attach coordinates for the convenience of output visualisation. Any coordinate information required by Reactions should be supplied as Variables. -```@meta -CurrentModule = PALEOboxes -``` -```@docs -NamedDimension -FixedCoord -``` - ## Variables ```@meta diff --git a/src/CoordsDims.jl b/src/CoordsDims.jl index 80bb662b..291efffc 100644 --- a/src/CoordsDims.jl +++ b/src/CoordsDims.jl @@ -37,12 +37,12 @@ function get_dimension end """ function set_coordinates!(obj, dimname, coordinates::Vector{String}) -Set coordinates attached to `dimname` for PALEO object `obj` +Set coordinates (variable names) attached to `dimname` for PALEO object `obj` -PALEO convention is that where possible `coordinates` contains three elements, for cell -midpoints, lower edges, upper edges, in that order. +PALEO convention is that where possible `coordinates` contains: +- three variable names, for cell midpoints, lower edges, upper edges, in that order. +- two variable names, for cell midpoints and a bounds variable (with a bounds dimension as the first dimensions) """ - function set_coordinates! end """ @@ -50,7 +50,8 @@ function set_coordinates! end Get coordinates (if any) attached to `dimname` for PALEO object `obj` -PALEO convention is that where possible `coordinates` contains three elements, for cell -midpoints, lower edges, upper edges, in that order. +PALEO convention is that where possible `coordinates` contains: +- three variable names, for cell midpoints, lower edges, upper edges, in that order. +- two variable names, for cell midpoints and a bounds variable (with a bounds dimension as the first dimensions) """ function get_coordinates end diff --git a/src/Grids.jl b/src/Grids.jl index 11f61d5f..dbf810b4 100644 --- a/src/Grids.jl +++ b/src/Grids.jl @@ -190,7 +190,14 @@ function get_subdomain(grid::PB.AbstractMesh, subdomainname::AbstractString) return subdomain end - + +""" + available_spaces(grid::PB.AbstractMesh) -> NTuple{PB.Space} + +Tuple of Spaces supported by this grid +""" +function available_spaces end + # generic handler when subdomainname present function PB.internal_size(space::Type{<:PB.AbstractSpace}, grid::PB.AbstractMeshOrNothing, subdomainname::AbstractString) @@ -428,15 +435,13 @@ Minimal Grid for a Vector Domain composed of columns (not necessarily forming a - `ncells::Int` total number of cells in this Domain - `Icolumns::Vector{Vector{Int}}`: Icolumns[n] should be the indices of column n, in order from surface to floor, where n is also the index of any associated boundary surface. -- `z_coords::Vector{FixedCoord}`: z coordinates of cell mid point, lower surface, upper surface -- `columnnames::Vector{Symbol}:` optional column names +- `columnnames::Vector{Symbol}`: optional column names +- `coordinates::Dict{String, Vector{String}}`: optional attached coordinates """ Base.@kwdef mutable struct UnstructuredColumnGrid <: PB.AbstractMesh ncells::Int64 Icolumns::Vector{Vector{Int}} - # z_coords::Vector{PB.FixedCoord} = PB.FixedCoord[] - "Define optional column names" columnnames::Vector{Symbol} = Symbol[] @@ -516,10 +521,8 @@ the PALEO internal representation (a Vector with a linear index) and a Cartesian # Fields - `ncells::Int64`: number of cells in Domain = `length(linear_index)` (may be subset of total points in `prod(dims)`) -- `dimnames::Vector{String}`: names of dimensions (ordered list) -- `dims::Vector{Int}`: sizes of dimensions (ordered list) -- `coords::Vector{Vector{Float64}}`: attached cell-centre coordinates for each dimension (ordered list) -- `coords_edges::Vector{Vector{Float64}}`: attached cell-edge coordinates for each dimension (ordered list) +- `dimensions::Vector{PB.NamedDimensions}`: names and sizes of cartesian spatial dimensions (ordered list) +- `coordinates::Dict{String, Vector{String}}`: optional attached coordinates """ Base.@kwdef mutable struct CartesianLinearGrid{N} <: PB.AbstractMesh @@ -591,11 +594,9 @@ _extra_dimensions(grid::CartesianLinearGrid) = grid.dimensions_extra nD grid with netcdf CF1.0 coordinates, using n-dimensional Arrays for PALEO internal representation of Variables # Fields -- `ncells::Int64`: number of cells in Domain = `length(linear_index)` (may be subset of total points in `prod(dims)`) -- `dimnames::Vector{String}`: names of dimensions (ordered list) -- `dims::Vector{Int}`: sizes of dimensions (ordered list) -- `coords::Vector{Vector{Float64}}`: attached cell-centre coordinates for each dimension (ordered list) -- `coords_edges::Vector{Vector{Float64}}`: attached cell-edge coordinates for each dimension (ordered list) +- `ncells::Int64`: number of cells in Domain = product of cartesian dimension sizes +- `dimensions::Vector{PB.NamedDimensions}`: names and sizes of cartesian spatial dimensions (ordered list) +- `coordinates::Dict{String, Vector{String}}`: optional attached coordinates """ Base.@kwdef mutable struct CartesianArrayGrid{N} <: PB.AbstractMesh From abb7ef6d32c2690fdbe260f687e2799346fac98a Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Mon, 30 Dec 2024 20:24:33 +0000 Subject: [PATCH 3/6] Small tidyups, remove CommonDataModel extension Move CommonDataModel extension changes to a separate PR --- Project cdm.toml | 71 ++++++ .../CommonDataModelExt/CommonDataModelExt.jl | 0 src/CoordsDims.jl | 6 +- src/Grids.jl | 52 ++-- src/Types cdm2.jl | 222 ++++++++++++++++++ src/Types.jl | 11 - 6 files changed, 336 insertions(+), 26 deletions(-) create mode 100644 Project cdm.toml rename {ext => ext_hide}/CommonDataModelExt/CommonDataModelExt.jl (100%) create mode 100644 src/Types cdm2.jl diff --git a/Project cdm.toml b/Project cdm.toml new file mode 100644 index 00000000..b3742098 --- /dev/null +++ b/Project cdm.toml @@ -0,0 +1,71 @@ +name = "PALEOboxes" +uuid = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4" +authors = ["Stuart Daines "] +version = "0.22.0" + +[deps] +Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MAT = "23992714-dd62-5051-b70f-ba57cb901cac" +NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +SIMD = "fdea26ae-647d-5447-a871-4b548cad5224" +SLEEF_jll = "63e82ce6-3d80-5af4-a84c-484b71ab98bb" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" + +[weakdeps] +CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" + +[extensions] +CommonDataModelExt = "CommonDataModel" + +[compat] +Atomix = "0.1, 1.0" +BenchmarkTools = "1.0" +CommonDataModel = "0.3.7" +DataFrames = "1.1" +DocStringExtensions = "0.8, 0.9" +Documenter = "1" +Graphs = "1.4" +Infiltrator = "1.0" +Interpolations = "0.13, 0.14, 0.15" +MAT = "0.10" +NCDatasets = "0.12, 0.13, 0.14" +OrderedCollections = "1.4" +PrecompileTools = "1.0" +Preferences = "1.2" +Revise = "3.1" +SIMD = "3.3" +SLEEF_jll = "3.4" +StaticArrays = "1.4" +StructArrays = "0.6, 0.7" +TestEnv = "1.0" +TimerOutputs = "0.5" +YAML = "0.4.7" +julia = "1.10" + +[extras] +CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["CommonDataModel", "Documenter", "Logging", "Test"] diff --git a/ext/CommonDataModelExt/CommonDataModelExt.jl b/ext_hide/CommonDataModelExt/CommonDataModelExt.jl similarity index 100% rename from ext/CommonDataModelExt/CommonDataModelExt.jl rename to ext_hide/CommonDataModelExt/CommonDataModelExt.jl diff --git a/src/CoordsDims.jl b/src/CoordsDims.jl index 291efffc..ee6a11e5 100644 --- a/src/CoordsDims.jl +++ b/src/CoordsDims.jl @@ -16,7 +16,11 @@ struct NamedDimension end function Base.show(io::IO, nd::NamedDimension) - print(io, "NamedDimension(name=", nd.name, ", size=", nd.size, ")") + if get(io, :typeinfo, nothing) === NamedDimension + print(io, "(name=", nd.name, ", size=", nd.size, ")") + else + print(io, "NamedDimension(name=", nd.name, ", size=", nd.size, ")") + end return nothing end diff --git a/src/Grids.jl b/src/Grids.jl index dbf810b4..99eaa5d7 100644 --- a/src/Grids.jl +++ b/src/Grids.jl @@ -684,24 +684,25 @@ function CartesianGrid( end """ - CartesianGrid(GridType, ncfilename::AbstractString, dimnames::Vector{<:AbstractString}; equalspacededges=false) - - > grid::GridType, coord_vars, coord_edges_vars + CartesianGrid(GridType, ncfilename::AbstractString, dimnames::Vector{<:AbstractString}; + equalspacededges=false, convert_edges_to_bounds=false) -> grid::GridType, coord_vars, coord_edges_bounds_vars Read a netcdf file with CF1.0 coordinates, and create corresponding a CartesianLinearGrid or CartesianArrayGrid from dimensions `dimnames`. -`coord_vars` and optional `coord_edges_vars` are Vectors of `coorddimname=>values` for creating coordinate variables. +`coord_vars` and optional `coord_edges_bounds_vars` are Vectors of `coorddimname => values` for creating coordinate variables. """ function CartesianGrid( GridType::Type{<:Union{CartesianLinearGrid, CartesianArrayGrid}}, ncfilename::AbstractString, dimnames::Vector{<:AbstractString}; - equalspacededges=false + equalspacededges=false, + convert_edges_to_bounds=false, ) @info "CartesianGrid creating $GridType{$(length(dimnames))} from dimnames=$dimnames in netcdf file $(ncfilename)" grid = GridType{length(dimnames)}() coord_vars = [] - coord_edges_vars = [] + coord_edges_bounds_vars = [] NCDatasets.Dataset(ncfilename) do ds # read dimensions and coordinates @@ -713,20 +714,43 @@ function CartesianGrid( grid.dimensions[i] = PB.NamedDimension(dimname, dim) v_cf = ds[dimname] coord_vals = Array(v_cf) + @assert ndims(coord_vals) == 1 + coord_vals = [cv for cv in coord_vals] # narrow type to eliminate unnecessary Missing push!(coord_vars, dimname=>coord_vals) if haskey(v_cf.attrib, "edges") edgesname = v_cf.attrib["edges"] - edgesdim = ds.dim[edgesname] + edgesdim = ds.dim[edgesname] @info " reading coordinate edges from $edgesname = $edgesdim" - push!(grid.dimensions_extra, PB.NamedDimension(edgesname, edgesdim)) - push!(coord_edges_vars, edgesname=>Array(ds[edgesname])) + edges_vals = Array(ds[edgesname]) + @assert size(edges_vals) == (edgesdim, ) + edges_vals = [ev for ev in edges_vals] + if convert_edges_to_bounds + boundsname = dimname*"_bnds" + @info " converting edges to bounds $boundsname" + coord_bounds = similar(coord_vals, eltype(coord_vals), 2, length(coord_vals)) + coord_bounds[1, :] = edges_vals[1:end-1] + coord_bounds[2, :] = edges_vals[2:end] + push!(coord_edges_bounds_vars, boundsname => coord_bounds) + else + push!(grid.dimensions_extra, PB.NamedDimension(edgesname, edgesdim)) + push!(coord_edges_bounds_vars, edgesname => edges_vals) + end elseif equalspacededges - edgesname = dimname*"_edges" ew = coord_vals[2] - coord_vals[1] - edgesvals = [coord_vals .- ew/2.0; coord_vals[end] + ew/2.0 ] - @info " assuming equal spacing $ew to calculate coordinate edges $edgesname" - push!(grid.dimensions_extra, PB.NamedDimension(edgesname, dim+1)) - push!(coord_edges_vars, edgesname=>edgesvals) + if convert_edges_to_bounds + boundsname = dimname*"_bnds" + @info " assuming equal spacing $ew to calculate coordinate bounds $boundsname" + coord_bounds = similar(coord_vals, eltype(coord_vals), 2, length(coord_vals)) + coord_bounds[1, :] = coord_vals .- ew/2.0 + coord_bounds[2, :] = coord_vals .+ ew/2.0 + push!(coord_edges_bounds_vars, boundsname => coord_bounds) + else + edgesname = dimname*"_edges" + @info " assuming equal spacing $ew to calculate coordinate edges $edgesname" + edges_vals = [coord_vals .- ew/2.0; coord_vals[end] + ew/2.0 ] + push!(grid.dimensions_extra, PB.NamedDimension(edgesname, dim + 1)) + push!(coord_edges_bounds_vars, edgesname => edges_vals) + end else @info " no edges attribute and equalspacededges=false" end @@ -760,7 +784,7 @@ function CartesianGrid( grid.ncells = prod(grid_size) end - return grid, coord_vars, coord_edges_vars + return grid, coord_vars, coord_edges_bounds_vars end """ diff --git a/src/Types cdm2.jl b/src/Types cdm2.jl new file mode 100644 index 00000000..94cc1c22 --- /dev/null +++ b/src/Types cdm2.jl @@ -0,0 +1,222 @@ +# TODO there doesn't seem to be an easy way of parameterising ModelData by an Array type ? +const PaleoArrayType = Array + +################################ +# CommonDataModel adaptor +################################ + +""" + CDModel(x) + +Create a CommonDataModel adaptor for PALEO object x +""" +function CDModel end + +################################ +# Parameters +############################### + +abstract type AbstractParameter +end + +################################# +# Reactions +################################### + + +""" + AbstractReactionMethod + +Base Type for Reaction Methods. + +""" +abstract type AbstractReactionMethod +end + + +abstract type AbstractReaction +end + + +############################################## +# Variables +############################################### + +""" + VariableBase + +A `Model` biogeochemical `Variable`. `Reaction`s access `Variable`s using derived Types [`VariableReaction`](@ref) +which are links to [`VariableDomain`](@ref)s. +""" +abstract type VariableBase +end + +""" + show_variables(obj, ...) -> Table + +Show all Variables attached to PALEO object `obj` +""" +function show_variables end + +""" + get_variable(obj, varname, ...) -> variable + +Get `variable <: VariableBase` by name from PALEO object +""" +function get_variable end + +""" + has_variable(obj, varname, ...) -> Bool + +True if PALEO object containts `varname` +""" +function has_variable end + + +""" + @enum VariableType + +Enumeration of `VariableBase` subtypes. Allowed values: +- `VariableReaction`: `VT_ReactProperty`, `VT_ReactDependency`, `VT_ReactContributor`, `VT_ReactTarget` +- `VariableDomain` : `VT_DomPropDep`, `VT_DomContribTarget` +""" +@enum VariableType::Cint begin + VT_Undefined = 0 + VT_ReactProperty + VT_ReactDependency + VT_ReactContributor + VT_ReactTarget + VT_DomPropDep + VT_DomContribTarget +end + + +abstract type VariableDomain <: VariableBase +end + +abstract type AbstractModel +end + +""" + AbstractDomain + +A model region containing Fields and Reactions that act on them. +""" +abstract type AbstractDomain +end + +abstract type AbstractSubdomain +end + +abstract type AbstractMesh +end + +const AbstractMeshOrNothing = Union{AbstractMesh, Nothing} + +""" + function get_mesh(obj, ...) + +Return an [`AbstractMesh`](@ref) for PALEO object `obj` +""" +function get_mesh end + +abstract type AbstractCellRange +end + +abstract type AbstractSpace +end + +abstract type AbstractData +end + +"parse eg \"ScalarData\" or \"PALEOboxes.ScalarData\" as ScalarData" +function Base.parse(::Type{AbstractData}, str::AbstractString) + dtype = tryparse(AbstractData, str) + + !isnothing(dtype) || + throw(ArgumentError("$str is not a subtype of AbstractData")) + return dtype +end + +function Base.tryparse(::Type{AbstractData}, str::AbstractString) + dtype = getproperty(@__MODULE__, Symbol(replace(str, "PALEOboxes."=>""))) + + return (dtype <: AbstractData) ? dtype : nothing +end + + +""" + AbstractField + +Defines a Field in a Domain +""" +abstract type AbstractField +end + +"`model` data arrays etc" +abstract type AbstractModelData +end + +"struct to hold Domain Field data" +abstract type AbstractDomainData +end + +""" + function get_table(obj, ...) + +Return a Tables.jl data table view for PALEO object `obj` +""" +function get_table end + + +############################################# +# ReactionMethodDispatchList +############################################# + +""" + ReactionMethodDispatchList + +Defines a list of [`ReactionMethod`](@ref) with corresponding [`CellRange`](@ref) +and views on Variable data (sub)arrays. +""" +struct ReactionMethodDispatchList{M <:Tuple, V <:Tuple, C <: Tuple} + methods::M + vardatas::V + cellranges::C +end + +ReactionMethodDispatchList(methods::Vector, vardatas::Vector, cellranges::Vector) = + ReactionMethodDispatchList(Tuple(methods), Tuple(vardatas), Tuple(cellranges)) + +struct ReactionMethodDispatchListNoGen + methods::Vector + vardatas::Vector + cellranges::Vector +end + +# See https://discourse.julialang.org/t/pretty-print-of-type/19555 +# Customize typeof function, as full type name is too verbose (as Tuples are of length ~ number of ReactionMethods to call) +Base.show(io::IO, ::Type{PALEOboxes.ReactionMethodDispatchList{M, V, C}}) where {M, V, C} = + print(io, "PALEOboxes.ReactionMethodDispatchList{M::Tuple, V::Tuple, C::Tuple each of length=$(fieldcount(M))}") + +function Base.show(io::IO, dispatchlist::ReactionMethodDispatchList) + print(io, typeof(dispatchlist)) +end + +function Base.show(io::IO, ::MIME"text/plain", dl::ReactionMethodDispatchList) + println(io, typeof(dl)) + for i in eachindex(dl.methodfns) + println(io, " ", dl.methods[i], ", ", dl.cellranges[i]) + end +end + +""" + infoerror(io::IOBuffer, message::AbstractString) + +Output accumulated log messages in io, then raise `ErrorException` with message +""" +function infoerror(io::IOBuffer, message::AbstractString) + s = String(take!(io)) + isempty(s) || @info s + error(message) +end \ No newline at end of file diff --git a/src/Types.jl b/src/Types.jl index 94cc1c22..8ed035ab 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -1,17 +1,6 @@ # TODO there doesn't seem to be an easy way of parameterising ModelData by an Array type ? const PaleoArrayType = Array -################################ -# CommonDataModel adaptor -################################ - -""" - CDModel(x) - -Create a CommonDataModel adaptor for PALEO object x -""" -function CDModel end - ################################ # Parameters ############################### From 2291de1cc7a33bc9cafe7ed3dca81d42206b1466 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Mon, 30 Dec 2024 20:33:39 +0000 Subject: [PATCH 4/6] really remove CommonDataModelExt --- Project cdm.toml | 71 ------ .../CommonDataModelExt/CommonDataModelExt.jl | 93 -------- src/Types cdm2.jl | 222 ------------------ 3 files changed, 386 deletions(-) delete mode 100644 Project cdm.toml delete mode 100644 ext_hide/CommonDataModelExt/CommonDataModelExt.jl delete mode 100644 src/Types cdm2.jl diff --git a/Project cdm.toml b/Project cdm.toml deleted file mode 100644 index b3742098..00000000 --- a/Project cdm.toml +++ /dev/null @@ -1,71 +0,0 @@ -name = "PALEOboxes" -uuid = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4" -authors = ["Stuart Daines "] -version = "0.22.0" - -[deps] -Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" -Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -MAT = "23992714-dd62-5051-b70f-ba57cb901cac" -NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" -SIMD = "fdea26ae-647d-5447-a871-4b548cad5224" -SLEEF_jll = "63e82ce6-3d80-5af4-a84c-484b71ab98bb" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" -TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" -TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" -YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" - -[weakdeps] -CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" - -[extensions] -CommonDataModelExt = "CommonDataModel" - -[compat] -Atomix = "0.1, 1.0" -BenchmarkTools = "1.0" -CommonDataModel = "0.3.7" -DataFrames = "1.1" -DocStringExtensions = "0.8, 0.9" -Documenter = "1" -Graphs = "1.4" -Infiltrator = "1.0" -Interpolations = "0.13, 0.14, 0.15" -MAT = "0.10" -NCDatasets = "0.12, 0.13, 0.14" -OrderedCollections = "1.4" -PrecompileTools = "1.0" -Preferences = "1.2" -Revise = "3.1" -SIMD = "3.3" -SLEEF_jll = "3.4" -StaticArrays = "1.4" -StructArrays = "0.6, 0.7" -TestEnv = "1.0" -TimerOutputs = "0.5" -YAML = "0.4.7" -julia = "1.10" - -[extras] -CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["CommonDataModel", "Documenter", "Logging", "Test"] diff --git a/ext_hide/CommonDataModelExt/CommonDataModelExt.jl b/ext_hide/CommonDataModelExt/CommonDataModelExt.jl deleted file mode 100644 index c2bdcfc4..00000000 --- a/ext_hide/CommonDataModelExt/CommonDataModelExt.jl +++ /dev/null @@ -1,93 +0,0 @@ -module CommonDataModelExt - -import PALEOboxes as PB -import CommonDataModel as CDM - -##################################### -# Wrapper types -###################################### - -struct ModelCDM <: CDM.AbstractDataset - model::PB.Model - modeldata::Union{Nothing, PB.AbstractModelData} -end - -PB.CDModel(model::PB.Model) = ModelCDM(model, nothing) -PB.CDModel(modeldata::PB.AbstractModelData) = ModelCDM(modeldata.model, modeldata) - -struct DomainCDM <: CDM.AbstractDataset - domain::PB.Domain - modeldata::Union{Nothing, PB.AbstractModelData} -end - -PB.CDModel(domain::PB.Domain) = DomainCDM(domain, nothing) - - -struct VariableDomainCDM <: CDM.AbstractDataset - variabledomain::PB.VariableDomain - modeldata::Union{Nothing, PB.AbstractModelData} - data -end - -PB.CDModel(variabledomain::PB.VariableDomain) = VariableDomainCDM(variabledomain, nothing, nothing) - -###################################### -# Model -###################################### - -# iterable with all group names -CDM.groupnames(m::ModelCDM) = [d.name for d in m.model.domains] - -CDM.group(m::ModelCDM, name::AbstractString) = DomainCDM(PB.get_domain(m.model, name; allow_not_found=false), m.modeldata) - -############################################### -# Domain -################################################ - -CDM.name(d::DomainCDM) = d.domain.name - -# TODO -# parentdataset(d::DomainCDM) - -# returns a list of variable names as strings -Base.keys(d::DomainCDM) = [v.name for v in PB.get_variables(d.domain)] - -function CDM.variable(d::DomainCDM, varname::AbstractString) - variabledomain = PB.get_variable(d.domain, varname; allow_not_found=false) - if isnothing(d.modeldata) - data = nothing - else - data = PB.get_data(variabledomain, d.modeldata) - end - return VariableDomainCDM(variabledomain, d.modeldata, data) -end - -CDM.dimnames(d::DomainCDM) = [nd.name for nd in PB.get_dimensions(d.domain)] - -CDM.dim(d::DomainCDM, name::AbstractString) = PB.get_dimension(d.domain, name).size - - -############################################### -# VariableDomain -################################################ - -CDM.name(v::VariableDomainCDM) = v.variabledomain.name - -CDM.dataset(v::VariableDomainCDM) = DomainCDM(v.variabledomain.domain, v.modeldata) - -CDM.dimnames(v::VariableDomainCDM) = [nd.name for nd in PB.get_dimensions(v.variabledomain)] - -Base.ndims(v::VariableDomainCDM) = length(CDM.dimnames(v)) - -Base.size(v::VariableDomainCDM) = (nd.size for nd in PB.get_dimensions(v.variabledomain)) - -CDM.attribnames(v::VariableDomainCDM) = keys(v.variabledomain.attributes) - -CDM.attrib(v::VariableDomainCDM, name::Symbol) = v.variabledomain.attributes[name] - -Base.getindex(v::VariableDomainCDM, indices...) = Base.getindex(v.data, indices...) - -Base.eltype(v::VariableDomainCDM) = Base.eltype(v.data) - - -end # module \ No newline at end of file diff --git a/src/Types cdm2.jl b/src/Types cdm2.jl deleted file mode 100644 index 94cc1c22..00000000 --- a/src/Types cdm2.jl +++ /dev/null @@ -1,222 +0,0 @@ -# TODO there doesn't seem to be an easy way of parameterising ModelData by an Array type ? -const PaleoArrayType = Array - -################################ -# CommonDataModel adaptor -################################ - -""" - CDModel(x) - -Create a CommonDataModel adaptor for PALEO object x -""" -function CDModel end - -################################ -# Parameters -############################### - -abstract type AbstractParameter -end - -################################# -# Reactions -################################### - - -""" - AbstractReactionMethod - -Base Type for Reaction Methods. - -""" -abstract type AbstractReactionMethod -end - - -abstract type AbstractReaction -end - - -############################################## -# Variables -############################################### - -""" - VariableBase - -A `Model` biogeochemical `Variable`. `Reaction`s access `Variable`s using derived Types [`VariableReaction`](@ref) -which are links to [`VariableDomain`](@ref)s. -""" -abstract type VariableBase -end - -""" - show_variables(obj, ...) -> Table - -Show all Variables attached to PALEO object `obj` -""" -function show_variables end - -""" - get_variable(obj, varname, ...) -> variable - -Get `variable <: VariableBase` by name from PALEO object -""" -function get_variable end - -""" - has_variable(obj, varname, ...) -> Bool - -True if PALEO object containts `varname` -""" -function has_variable end - - -""" - @enum VariableType - -Enumeration of `VariableBase` subtypes. Allowed values: -- `VariableReaction`: `VT_ReactProperty`, `VT_ReactDependency`, `VT_ReactContributor`, `VT_ReactTarget` -- `VariableDomain` : `VT_DomPropDep`, `VT_DomContribTarget` -""" -@enum VariableType::Cint begin - VT_Undefined = 0 - VT_ReactProperty - VT_ReactDependency - VT_ReactContributor - VT_ReactTarget - VT_DomPropDep - VT_DomContribTarget -end - - -abstract type VariableDomain <: VariableBase -end - -abstract type AbstractModel -end - -""" - AbstractDomain - -A model region containing Fields and Reactions that act on them. -""" -abstract type AbstractDomain -end - -abstract type AbstractSubdomain -end - -abstract type AbstractMesh -end - -const AbstractMeshOrNothing = Union{AbstractMesh, Nothing} - -""" - function get_mesh(obj, ...) - -Return an [`AbstractMesh`](@ref) for PALEO object `obj` -""" -function get_mesh end - -abstract type AbstractCellRange -end - -abstract type AbstractSpace -end - -abstract type AbstractData -end - -"parse eg \"ScalarData\" or \"PALEOboxes.ScalarData\" as ScalarData" -function Base.parse(::Type{AbstractData}, str::AbstractString) - dtype = tryparse(AbstractData, str) - - !isnothing(dtype) || - throw(ArgumentError("$str is not a subtype of AbstractData")) - return dtype -end - -function Base.tryparse(::Type{AbstractData}, str::AbstractString) - dtype = getproperty(@__MODULE__, Symbol(replace(str, "PALEOboxes."=>""))) - - return (dtype <: AbstractData) ? dtype : nothing -end - - -""" - AbstractField - -Defines a Field in a Domain -""" -abstract type AbstractField -end - -"`model` data arrays etc" -abstract type AbstractModelData -end - -"struct to hold Domain Field data" -abstract type AbstractDomainData -end - -""" - function get_table(obj, ...) - -Return a Tables.jl data table view for PALEO object `obj` -""" -function get_table end - - -############################################# -# ReactionMethodDispatchList -############################################# - -""" - ReactionMethodDispatchList - -Defines a list of [`ReactionMethod`](@ref) with corresponding [`CellRange`](@ref) -and views on Variable data (sub)arrays. -""" -struct ReactionMethodDispatchList{M <:Tuple, V <:Tuple, C <: Tuple} - methods::M - vardatas::V - cellranges::C -end - -ReactionMethodDispatchList(methods::Vector, vardatas::Vector, cellranges::Vector) = - ReactionMethodDispatchList(Tuple(methods), Tuple(vardatas), Tuple(cellranges)) - -struct ReactionMethodDispatchListNoGen - methods::Vector - vardatas::Vector - cellranges::Vector -end - -# See https://discourse.julialang.org/t/pretty-print-of-type/19555 -# Customize typeof function, as full type name is too verbose (as Tuples are of length ~ number of ReactionMethods to call) -Base.show(io::IO, ::Type{PALEOboxes.ReactionMethodDispatchList{M, V, C}}) where {M, V, C} = - print(io, "PALEOboxes.ReactionMethodDispatchList{M::Tuple, V::Tuple, C::Tuple each of length=$(fieldcount(M))}") - -function Base.show(io::IO, dispatchlist::ReactionMethodDispatchList) - print(io, typeof(dispatchlist)) -end - -function Base.show(io::IO, ::MIME"text/plain", dl::ReactionMethodDispatchList) - println(io, typeof(dl)) - for i in eachindex(dl.methodfns) - println(io, " ", dl.methods[i], ", ", dl.cellranges[i]) - end -end - -""" - infoerror(io::IOBuffer, message::AbstractString) - -Output accumulated log messages in io, then raise `ErrorException` with message -""" -function infoerror(io::IOBuffer, message::AbstractString) - s = String(take!(io)) - isempty(s) || @info s - error(message) -end \ No newline at end of file From ed1c6a0cdc72bc2b6fad9700da3249102e5a0ce8 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Mon, 30 Dec 2024 20:35:16 +0000 Subject: [PATCH 5/6] really remove CommonDataModelExt --- Project.toml | 7 ------- 1 file changed, 7 deletions(-) diff --git a/Project.toml b/Project.toml index b3742098..c8aa7d4e 100644 --- a/Project.toml +++ b/Project.toml @@ -30,16 +30,9 @@ TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" -[weakdeps] -CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" - -[extensions] -CommonDataModelExt = "CommonDataModel" - [compat] Atomix = "0.1, 1.0" BenchmarkTools = "1.0" -CommonDataModel = "0.3.7" DataFrames = "1.1" DocStringExtensions = "0.8, 0.9" Documenter = "1" From 6cb1c14b50299d4bf4fb0b7360e154077f33115e Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Mon, 30 Dec 2024 20:38:01 +0000 Subject: [PATCH 6/6] fix Project.toml --- Project.toml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index c8aa7d4e..b4157b05 100644 --- a/Project.toml +++ b/Project.toml @@ -52,13 +52,12 @@ StructArrays = "0.6, 0.7" TestEnv = "1.0" TimerOutputs = "0.5" YAML = "0.4.7" -julia = "1.10" +julia = "1.6" [extras] -CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["CommonDataModel", "Documenter", "Logging", "Test"] +test = ["Documenter", "Logging", "Test"]