From 8e859f005f172db94b53f624f2cb2d5baf37f1da Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Sun, 2 Jul 2023 16:32:04 -0700 Subject: [PATCH 1/2] Separate model components into new files --- src/simulation_components/atmos.jl | 66 +++ .../domains.jl} | 553 +----------------- src/{ => simulation_components}/floe.jl | 0 src/simulation_components/grids.jl | 206 +++++++ src/simulation_components/model.jl | 117 ++++ src/simulation_components/oceans.jl | 171 ++++++ src/{ => simulation_components}/simulation.jl | 0 7 files changed, 561 insertions(+), 552 deletions(-) create mode 100644 src/simulation_components/atmos.jl rename src/{model.jl => simulation_components/domains.jl} (54%) rename src/{ => simulation_components}/floe.jl (100%) create mode 100644 src/simulation_components/grids.jl create mode 100644 src/simulation_components/model.jl create mode 100644 src/simulation_components/oceans.jl rename src/{ => simulation_components}/simulation.jl (100%) diff --git a/src/simulation_components/atmos.jl b/src/simulation_components/atmos.jl new file mode 100644 index 0000000..b1f0d7b --- /dev/null +++ b/src/simulation_components/atmos.jl @@ -0,0 +1,66 @@ +""" +Structs and functions used to define a Subzero atmosphere +""" + + +""" +Atmos velocities in the x-direction (u) and y-direction (v). u and v should +match the size of the corresponding model grid so that there is one x and y +velocity value for each grid cell. Atmos also needs temperature at the +atmosphere/ice interface in each grid cell. Model cannot be constructed if size +of atmos fields and grid do not match. +""" +struct Atmos{FT<:AbstractFloat} + u::Matrix{FT} + v::Matrix{FT} + temp::Matrix{FT} + function Atmos{FT}(u, v, temp) where {FT <: AbstractFloat} + if !(size(u) == size(v)) + throw(ArgumentError("One or more of the atmosphere vector fields \ + aren't the same dimension.")) + end + new{FT}(u, v, temp) + end +end + +""" + Atmos(::Type{FT}, args...) + +A float type FT can be provided as the first argument of any Atmos constructor. +An Atmos of type FT will be created by passing all other arguments to the +correct constructor. +""" +Atmos(::Type{FT}, args...) where {FT <: AbstractFloat} = + Atmos{FT}(args...) + +""" + Atmos(args...) + +If a type isn't specified, Atmos will be of type Float64 and the correct +constructor will be called with all other arguments. +""" +Atmos(args...) = Atmos{Float64}(args...) + +""" + Atmos{FT}(grid, u, v) + +Construct model atmosphere of type FT. +Inputs: + grid model's grid + u Atmos x-velocity for each grid cell + v Atmos y-velocity for each grid cell + temp temperature at atmopshere/ice interface per grid cell +Output: + Atmosphere of type FT with constant velocity and temperature over domain. +""" +Atmos{FT}( + grid::AbstractGrid, + u, + v, + temp, +) where {FT <: AbstractFloat} = + return Atmos{FT}( + fill(FT(u), grid.Nx + 1, grid.Ny + 1), + fill(FT(v), grid.Nx + 1, grid.Ny + 1), + fill(FT(temp), grid.Nx + 1, grid.Ny + 1), + ) diff --git a/src/model.jl b/src/simulation_components/domains.jl similarity index 54% rename from src/model.jl rename to src/simulation_components/domains.jl index 77e67c0..0ca64d9 100644 --- a/src/model.jl +++ b/src/simulation_components/domains.jl @@ -1,443 +1,7 @@ """ -Structs and functions used to define a Subzero model +Structs and functions used to define a Subzero domains """ -### ---------------- Grid Fields, Structs, and Constructors ---------------- ### -""" - CellFloes{FT<:AbstractFloat} - -Struct that tracks which floes are within given cell, as well as the translation -vector needed to move floe each from current postion into cell if it is in cell -due to periodic boundaries. Each index in floeidx is the index of a floe within -the cell and the Δx and Δy with the same index are that floe's translation -vector. -Note: the floeidx is the index of grid cells centered on grid lines, not on the -grid cells defined by the regular, rectilinear grid. -""" -struct CellFloes{FT<:AbstractFloat} - floeidx::Vector{Int} - Δx::Vector{FT} - Δy::Vector{FT} -end - -""" - CellFloes{FT}() - -Constructs an CellFloes object with empty lists for fields. -""" -CellFloes{FT}() where {FT} = CellFloes{FT}( - Vector{Int}(), - Vector{FT}(), - Vector{FT}() -) -""" - empty!(cell::CellFloes) - -Empties the vectors within a CellFloes object -""" -function Base.empty!(cell::CellFloes) - empty!(cell.floeidx) - empty!(cell.Δx) - empty!(cell.Δy) - return -end - -""" - AbstractGrid - -An abstract type for the grid that model will be simulated on. -Affects calculation on the grid. -""" -abstract type AbstractGrid{FT<:AbstractFloat} end - -""" - RegRectilinearGrid{FT<:AbstractFloat}<:AbstractGrid - -Tessellation of 2-dimensional Euclidean space into n-by-m congruent rectangles. -- Nx: number of grid cells in the x-direction -- Ny: number of grid cells in the y-direction -- x0: value of first x grid line -- xf: value of final x grid line -- y0: value of first y grid line -- yf: value of final y grid line -- Δx: grid cell width -- Δy: grid cell height -- floe_locations: matrix of CellFloes, which keep track of which floes are in - each cell -""" -struct RegRectilinearGrid{FT}<:AbstractGrid{FT} - Nx::Int - Ny::Int - x0::FT - xf::FT - y0::FT - yf::FT - Δx::FT - Δy::FT - floe_locations::Matrix{CellFloes{FT}} - - function RegRectilinearGrid{FT}( - Nx, - Ny, - x0, - xf, - y0, - yf, - Δx, - Δy, - floe_locations, - ) where {FT} - if size(floe_locations) != (Nx + 1, Ny + 1) - throw(ArgumentError("Floe location matrix needs to be the same \ - dimensions as grid lines.")) - end - if (xf - x0) / Δx != Nx - throw(ArgumentError("X-grid extent and grid cell width don't match \ - with number of grid cells in x-direction.")) - end - if (yf - y0) / Δy != Ny - throw(ArgumentError("Y-grid extent and grid cell height don't match \ - with number of grid cells in y-direction.")) - end - return new{FT}(Nx, Ny, x0, xf, y0, yf, Δx, Δy, floe_locations) - end -end - -""" - RegRectilinearGrid(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any RegRectilinearGrid -constructor. A RegRectilinearGrid of type FT will be created by passing all -other arguments to the correct constructor. -""" -RegRectilinearGrid(::Type{FT}, args...) where {FT <: AbstractFloat} = - RegRectilinearGrid{FT}(args...) - -""" - RegRectilinearGrid(args...) - -If a type isn't specified, RegRectilinearGrid will be of type Float64 and the -correct constructor will be called with all other arguments. -""" -RegRectilinearGrid(args...) = RegRectilinearGrid{Float64}(args...) - -""" - RegRectilinearGrid{FT}( - xbounds::Tuple, - ybounds::Tuple, - Δx, - Δy, - ) - -Construct a RegRectilinearGrid for model given bounds for grid x and y and grid -cell dimensions in meters. -Inputs: - xbounds bound of grid x-direction in form (left, right) - ybounds bound of grid y-direction in form (bottom, top) - Δx length/height of grid cells in x-direction - Δy length/height of grid cells in y-direction -Output: - RegRectilinearGrid from lx to ux and height from ly to uy with grid squares - of size Δx by Δy -Warning: - If Δx doesn't evenly divide x length (lu-lx) or Δy doesn't evenly divide y - length (uy-ly) you won't get full size grid. The grid will be "trimmed" to - the nearest full grid square in both directions. -""" -function RegRectilinearGrid{FT}( - xbounds::Tuple, - ybounds::Tuple, - Δx, - Δy, -) where {FT <: AbstractFloat} - Nx = floor(Int, (xbounds[2] - xbounds[1]) / Δx) - Ny = floor(Int, (ybounds[2] - ybounds[1]) / Δy) - return RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds[1], - xbounds[1] + Nx * Δx, - ybounds[1], - ybounds[1] + Ny * Δy, - Δx, - Δy, - [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], - ) -end - -get_grid_lines(g0, gf, Δg) = g0:Δg:gf -get_grid_centers(g0, gf, Δg) = (g0 + Δg/2):Δg:(gf - Δg/2) - -""" - RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds::Tuple{Real, Real}, - ybounds::Tuple{Real, Real}, - ) where {FT <: AbstractFloat} - -Construct a RegRectilinearGrid for model given bounds for grid x and y and the -number of grid cells in both the x and y direction. -Inputs: - Nx number of grid cells in the x-direction - Ny number of grid cells in the y-direction - xbounds bound of grid x-direction in form (left, right) - ybounds bound of grid y-direction in form (bottom, top) -Output: - RegRectilinearGrid with width and height determined by xbound and ybounds - and the number of grid cells in the x-direction and y-direction determined - by dims. -""" -RegRectilinearGrid{FT}( - Nx::Int, - Ny::Int, - xbounds::Tuple, - ybounds::Tuple, -) where {FT <: AbstractFloat} = - RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds[1], - xbounds[2], - ybounds[1], - ybounds[2], - (xbounds[2] - xbounds[1]) / Nx, - (ybounds[2] - ybounds[1]) / Ny, - [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], - ) - -### --------------- Ocean Fields, Structs, and Constructors --------------- ### -""" - IceStressCell{FT<:AbstractFloat} - -Struct to collect stress from ice floes on ocean grid cells. One IceStressCell -corresponds to one grid cell. It holds a list of running totals of stress on the -cell, and a running list of the number of points making up those running totals. -Each element in the list corresponds to one floe, which is denoted in the -corresponding CellFloes matrix within the grid. -""" -struct IceStressCell{FT<:AbstractFloat} - τx::Vector{FT} - τy::Vector{FT} - npoints::Vector{Int} -end - -""" - IceStressCell{FT}() - -Constructs an IceStressCell object with empty lists for fields. -""" -IceStressCell{FT}() where {FT} = IceStressCell{FT}( - Vector{FT}(), - Vector{FT}(), - Vector{Int}() -) -""" - empty!(cell::IceStressCell) - -Empties the vectors within an IceStressCell -""" -function Base.empty!(cell::IceStressCell) - empty!(cell.τx) - empty!(cell.τy) - empty!(cell.npoints) - return -end - -""" - Ocean{FT<:AbstractFloat} - -Simulation ocean holding ocean values on the grid-scale with matricies of the -same size as the model's grid. The struct has the following fields: -- u is the ocean velocities in the x-direction for each grid cell -- v is the ocean velocities in the y-direction for each grid cell -- temp is the ocean temperature for each grid cell -- hflx_factor is a factor to calculate the ocean-atmosphere heat flux for a - cell in that grid cell by multiplying by its height -- si_frac is the fraction of area in each grid cell that is covered in sea-ice - -Ocean fields must all be matricies with dimensions equal to the number of grid -lines in the model's grid. -Note: If a periodic boundary is used in the domain, the last grid cell in that -direction will not be used as it is equivalent to the first grid cell. Thus, for -all of these fields, the first and last value in the x and/or y direction should -be equal if the east-west or north-south boundary pair are periodic -respectively. -""" -struct Ocean{FT<:AbstractFloat} - u::Matrix{FT} - v::Matrix{FT} - temp::Matrix{FT} - hflx_factor::Matrix{FT} - scells::Matrix{IceStressCell{FT}} - τx::Matrix{FT} - τy::Matrix{FT} - si_frac::Matrix{FT} - dissolved::Matrix{FT} - - function Ocean{FT}( - u, - v, - temp, - hflx, - scells, - τx, - τy, - si_frac, - dissolved, - ) where {FT <: AbstractFloat} - if !all(-3 .<= temp .<= 0) - @warn "Ocean temperatures are above the range for freezing. The \ - thermodynamics aren't currently setup for these conditions." - end - if !(size(u) == size(v) == size(τx) == size(τy)) - throw(ArgumentError("One or more of the ocean vector fields aren't \ - the same dimension.")) - end - if !(size(temp) == size(hflx) == size(si_frac) == size(dissolved)) - throw(ArgumentError("One or more of the ocean tracer fields aren't \ - the same dimension.")) - end - new{FT}(u, v, temp, hflx, scells, τx, τy, si_frac, dissolved) - end -end -""" - Ocean(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any Ocean constructor. -An Ocean of type FT will be created by passing all other arguments to the -correct constructor. -""" -Ocean(::Type{FT}, args...) where {FT <: AbstractFloat} = - Ocean{FT}(args...) - -""" - Ocean(args...) - -If a type isn't specified, Ocean will be of type Float64 and the correct -constructor will be called with all other arguments. -""" -Ocean(args...) = Ocean{Float64}(args...) - -""" - Ocean{FT}(u, v, temp) - -Construct model ocean. -Inputs: - u ocean x-velocity matrix with u for each grid line - v ocean y-velocity matrix with u for each grid line - temp temperature matrix with ocean/ice interface temperature for - each grid line -Output: - Model ocean with given velocity and temperature fields on each grid line. -""" -function Ocean{FT}( - u, - v, - temp, -) where {FT <: AbstractFloat} - Nx, Ny = size(temp) .- 1 - return Ocean{FT}( - u, - v, - temp, - zeros(FT, Nx + 1, Ny + 1), # heat flux - [IceStressCell{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], - zeros(FT, Nx + 1, Ny + 1), # x-stress - zeros(FT, Nx + 1, Ny + 1), # y-stress - zeros(FT, Nx + 1, Ny + 1), # sea ice fraction - zeros(FT, Nx + 1, Ny + 1), # dissolved - ) -end - -""" - Ocean{FT}(grid, u, v, temp) - -Construct model's ocean. -Inputs: - grid model grid - u ocean x-velocity for each grid line - v ocean y-velocity for each grid line - temp temperature at ocean/ice interface per grid cell -Output: - Ocean with constant velocity and temperature on each grid line. -""" -Ocean{FT}( - grid::AbstractGrid, - u, - v, - temp, -) where {FT <: AbstractFloat} = - Ocean{FT}( - fill(FT(u), grid.Nx + 1, grid.Ny + 1), - fill(FT(v), grid.Nx + 1, grid.Ny + 1), - fill(FT(temp), grid.Nx + 1, grid.Ny + 1), - ) - -### ------------- Atmosphere Fields, Structs, and Constructors ------------- ### -""" -Atmos velocities in the x-direction (u) and y-direction (v). u and v should -match the size of the corresponding model grid so that there is one x and y -velocity value for each grid cell. Atmos also needs temperature at the -atmosphere/ice interface in each grid cell. Model cannot be constructed if size -of atmos fields and grid do not match. -""" -struct Atmos{FT<:AbstractFloat} - u::Matrix{FT} - v::Matrix{FT} - temp::Matrix{FT} - function Atmos{FT}(u, v, temp) where {FT <: AbstractFloat} - if !(size(u) == size(v)) - throw(ArgumentError("One or more of the atmosphere vector fields \ - aren't the same dimension.")) - end - new{FT}(u, v, temp) - end -end - -""" - Atmos(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any Atmos constructor. -An Atmos of type FT will be created by passing all other arguments to the -correct constructor. -""" -Atmos(::Type{FT}, args...) where {FT <: AbstractFloat} = - Atmos{FT}(args...) - -""" - Atmos(args...) - -If a type isn't specified, Atmos will be of type Float64 and the correct -constructor will be called with all other arguments. -""" -Atmos(args...) = Atmos{Float64}(args...) - -""" - Atmos{FT}(grid, u, v) - -Construct model atmosphere of type FT. -Inputs: - grid model's grid - u Atmos x-velocity for each grid cell - v Atmos y-velocity for each grid cell - temp temperature at atmopshere/ice interface per grid cell -Output: - Atmosphere of type FT with constant velocity and temperature over domain. -""" -Atmos{FT}( - grid::AbstractGrid, - u, - v, - temp, -) where {FT <: AbstractFloat} = - return Atmos{FT}( - fill(FT(u), grid.Nx + 1, grid.Ny + 1), - fill(FT(v), grid.Nx + 1, grid.Ny + 1), - fill(FT(temp), grid.Nx + 1, grid.Ny + 1), - ) - -### --------------- Domain Fields, Structs, and Constructors --------------- ### """ AbstractDirection @@ -1102,118 +666,3 @@ Domain( west, StructArray{TopographyElement{FT}}(undef, 0), ) - -### --------------- Model Fields, Structs, and Constructors --------------- ### -""" - domain_in_grid(domain, grid) - -Checks if given rectangular domain is within given grid and gives user a warning -if domain is not of maximum possible size given grid dimensions. -Inputs: - domain - grid -Outputs: - true if domain is within grid bounds, else false -""" -function domain_in_grid(domain::Domain, grid::AbstractGrid) - northval = domain.north.val - southval = domain.south.val - eastval = domain.east.val - westval = domain.west.val - if (northval <= grid.yf && - southval >= grid.y0 && - eastval <= grid.xf && - westval >= grid.x0) - if (northval != grid.yf || - southval != grid.y0 || - eastval != grid.xf || - westval != grid.x0) - @warn "At least one wall of domain is smaller than grid. This \ - could lead to unneeded computation. Consider making grid \ - smaller or domain larger." - end - return true - end - return false -end - -""" -Model which holds grid, ocean, atmos structs, each with the same underlying -float type (either Float32 of Float64) and size. It also holds the domain -information, which includes the topography and the boundaries. It holds a -StructArray of floe structs, again each relying on the same underlying float -type. Finally, it also holds the maximum floe id used thus far in the -simulation. This should be the length of the floes array at the beginning of the -run. -""" -struct Model{ - FT<:AbstractFloat, - GT<:AbstractGrid{FT}, - DT<:Domain{ - FT, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - }, -} - grid::GT - ocean::Ocean{FT} - atmos::Atmos{FT} - domain::DT - floes::StructArray{Floe{FT}} # See floes.jl for floe creation - - function Model{FT, GT, DT}( - grid::GT, - ocean::Ocean{FT}, - atmos::Atmos{FT}, - domain::DT, - floes::StructArray{Floe{FT}}, - ) where { - FT<:AbstractFloat, - GT<:AbstractGrid{FT}, - DT<:Domain{ - FT, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - }, - } - if !domain_in_grid(domain, grid) - throw(ArgumentError("Domain does not fit within grid.")) - elseif ( - size(ocean.u) != size(atmos.u) || - size(ocean.v) != size(atmos.v) || - size(ocean.temp) != size(atmos.temp) - ) - throw(ArgumentError("Ocean and atmosphere are not on the same grid.\ - This is not supported yet.")) - end - if any(ocean.temp .< atmos.temp) - @warn "In at least one grid cell the atmosphere temperature is \ - warmer than the ocean. This is not a situation in which the \ - thermodynamics are setup for right now." - end - new{FT, GT, DT}(grid, ocean, atmos, domain, floes) - end - - Model( - grid::GT, - ocean::Ocean{FT}, - atmos::Atmos{FT}, - domain::DT, - floes::StructArray{Floe{FT}}, - ) where { - FT<:AbstractFloat, - GT<:AbstractGrid{FT}, - DT<:Domain{ - FT, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - <:AbstractBoundary, - }, - } = - Model{FT, GT, DT}(grid, ocean, atmos, domain, floes) -end \ No newline at end of file diff --git a/src/floe.jl b/src/simulation_components/floe.jl similarity index 100% rename from src/floe.jl rename to src/simulation_components/floe.jl diff --git a/src/simulation_components/grids.jl b/src/simulation_components/grids.jl new file mode 100644 index 0000000..2a60ab7 --- /dev/null +++ b/src/simulation_components/grids.jl @@ -0,0 +1,206 @@ +""" +Structs and functions used to define a Subzero grids +""" + +""" + CellFloes{FT<:AbstractFloat} + +Struct that tracks which floes are within given cell, as well as the translation +vector needed to move floe each from current postion into cell if it is in cell +due to periodic boundaries. Each index in floeidx is the index of a floe within +the cell and the Δx and Δy with the same index are that floe's translation +vector. +Note: the floeidx is the index of grid cells centered on grid lines, not on the +grid cells defined by the regular, rectilinear grid. +""" +struct CellFloes{FT<:AbstractFloat} + floeidx::Vector{Int} + Δx::Vector{FT} + Δy::Vector{FT} +end + +""" + CellFloes{FT}() + +Constructs an CellFloes object with empty lists for fields. +""" +CellFloes{FT}() where {FT} = CellFloes{FT}( + Vector{Int}(), + Vector{FT}(), + Vector{FT}() +) +""" + empty!(cell::CellFloes) + +Empties the vectors within a CellFloes object +""" +function Base.empty!(cell::CellFloes) + empty!(cell.floeidx) + empty!(cell.Δx) + empty!(cell.Δy) + return +end + +""" + AbstractGrid + +An abstract type for the grid that model will be simulated on. +Affects calculation on the grid. +""" +abstract type AbstractGrid{FT<:AbstractFloat} end + +""" + RegRectilinearGrid{FT<:AbstractFloat}<:AbstractGrid + +Tessellation of 2-dimensional Euclidean space into n-by-m congruent rectangles. +- Nx: number of grid cells in the x-direction +- Ny: number of grid cells in the y-direction +- x0: value of first x grid line +- xf: value of final x grid line +- y0: value of first y grid line +- yf: value of final y grid line +- Δx: grid cell width +- Δy: grid cell height +- floe_locations: matrix of CellFloes, which keep track of which floes are in + each cell +""" +struct RegRectilinearGrid{FT}<:AbstractGrid{FT} + Nx::Int + Ny::Int + x0::FT + xf::FT + y0::FT + yf::FT + Δx::FT + Δy::FT + floe_locations::Matrix{CellFloes{FT}} + + function RegRectilinearGrid{FT}( + Nx, + Ny, + x0, + xf, + y0, + yf, + Δx, + Δy, + floe_locations, + ) where {FT} + if size(floe_locations) != (Nx + 1, Ny + 1) + throw(ArgumentError("Floe location matrix needs to be the same \ + dimensions as grid lines.")) + end + if (xf - x0) / Δx != Nx + throw(ArgumentError("X-grid extent and grid cell width don't match \ + with number of grid cells in x-direction.")) + end + if (yf - y0) / Δy != Ny + throw(ArgumentError("Y-grid extent and grid cell height don't match \ + with number of grid cells in y-direction.")) + end + return new{FT}(Nx, Ny, x0, xf, y0, yf, Δx, Δy, floe_locations) + end +end + +""" + RegRectilinearGrid(::Type{FT}, args...) + +A float type FT can be provided as the first argument of any RegRectilinearGrid +constructor. A RegRectilinearGrid of type FT will be created by passing all +other arguments to the correct constructor. +""" +RegRectilinearGrid(::Type{FT}, args...) where {FT <: AbstractFloat} = + RegRectilinearGrid{FT}(args...) + +""" + RegRectilinearGrid(args...) + +If a type isn't specified, RegRectilinearGrid will be of type Float64 and the +correct constructor will be called with all other arguments. +""" +RegRectilinearGrid(args...) = RegRectilinearGrid{Float64}(args...) + +""" + RegRectilinearGrid{FT}( + xbounds::Tuple, + ybounds::Tuple, + Δx, + Δy, + ) + +Construct a RegRectilinearGrid for model given bounds for grid x and y and grid +cell dimensions in meters. +Inputs: + xbounds bound of grid x-direction in form (left, right) + ybounds bound of grid y-direction in form (bottom, top) + Δx length/height of grid cells in x-direction + Δy length/height of grid cells in y-direction +Output: + RegRectilinearGrid from lx to ux and height from ly to uy with grid squares + of size Δx by Δy +Warning: + If Δx doesn't evenly divide x length (lu-lx) or Δy doesn't evenly divide y + length (uy-ly) you won't get full size grid. The grid will be "trimmed" to + the nearest full grid square in both directions. +""" +function RegRectilinearGrid{FT}( + xbounds::Tuple, + ybounds::Tuple, + Δx, + Δy, +) where {FT <: AbstractFloat} + Nx = floor(Int, (xbounds[2] - xbounds[1]) / Δx) + Ny = floor(Int, (ybounds[2] - ybounds[1]) / Δy) + return RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds[1], + xbounds[1] + Nx * Δx, + ybounds[1], + ybounds[1] + Ny * Δy, + Δx, + Δy, + [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], + ) +end + +get_grid_lines(g0, gf, Δg) = g0:Δg:gf +get_grid_centers(g0, gf, Δg) = (g0 + Δg/2):Δg:(gf - Δg/2) + +""" + RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds::Tuple{Real, Real}, + ybounds::Tuple{Real, Real}, + ) where {FT <: AbstractFloat} + +Construct a RegRectilinearGrid for model given bounds for grid x and y and the +number of grid cells in both the x and y direction. +Inputs: + Nx number of grid cells in the x-direction + Ny number of grid cells in the y-direction + xbounds bound of grid x-direction in form (left, right) + ybounds bound of grid y-direction in form (bottom, top) +Output: + RegRectilinearGrid with width and height determined by xbound and ybounds + and the number of grid cells in the x-direction and y-direction determined + by dims. +""" +RegRectilinearGrid{FT}( + Nx::Int, + Ny::Int, + xbounds::Tuple, + ybounds::Tuple, +) where {FT <: AbstractFloat} = + RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds[1], + xbounds[2], + ybounds[1], + ybounds[2], + (xbounds[2] - xbounds[1]) / Nx, + (ybounds[2] - ybounds[1]) / Ny, + [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], + ) diff --git a/src/simulation_components/model.jl b/src/simulation_components/model.jl new file mode 100644 index 0000000..0a9e8a0 --- /dev/null +++ b/src/simulation_components/model.jl @@ -0,0 +1,117 @@ +""" +Structs and functions used to define a Subzero model +""" + +""" + domain_in_grid(domain, grid) + +Checks if given rectangular domain is within given grid and gives user a warning +if domain is not of maximum possible size given grid dimensions. +Inputs: + domain + grid +Outputs: + true if domain is within grid bounds, else false +""" +function domain_in_grid(domain::Domain, grid::AbstractGrid) + northval = domain.north.val + southval = domain.south.val + eastval = domain.east.val + westval = domain.west.val + if (northval <= grid.yf && + southval >= grid.y0 && + eastval <= grid.xf && + westval >= grid.x0) + if (northval != grid.yf || + southval != grid.y0 || + eastval != grid.xf || + westval != grid.x0) + @warn "At least one wall of domain is smaller than grid. This \ + could lead to unneeded computation. Consider making grid \ + smaller or domain larger." + end + return true + end + return false +end + +""" +Model which holds grid, ocean, atmos structs, each with the same underlying +float type (either Float32 of Float64) and size. It also holds the domain +information, which includes the topography and the boundaries. It holds a +StructArray of floe structs, again each relying on the same underlying float +type. Finally, it also holds the maximum floe id used thus far in the +simulation. This should be the length of the floes array at the beginning of the +run. +""" +struct Model{ + FT<:AbstractFloat, + GT<:AbstractGrid{FT}, + DT<:Domain{ + FT, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + }, +} + grid::GT + ocean::Ocean{FT} + atmos::Atmos{FT} + domain::DT + floes::StructArray{Floe{FT}} # See floes.jl for floe creation + + function Model{FT, GT, DT}( + grid::GT, + ocean::Ocean{FT}, + atmos::Atmos{FT}, + domain::DT, + floes::StructArray{Floe{FT}}, + ) where { + FT<:AbstractFloat, + GT<:AbstractGrid{FT}, + DT<:Domain{ + FT, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + }, + } + if !domain_in_grid(domain, grid) + throw(ArgumentError("Domain does not fit within grid.")) + elseif ( + size(ocean.u) != size(atmos.u) || + size(ocean.v) != size(atmos.v) || + size(ocean.temp) != size(atmos.temp) + ) + throw(ArgumentError("Ocean and atmosphere are not on the same grid.\ + This is not supported yet.")) + end + if any(ocean.temp .< atmos.temp) + @warn "In at least one grid cell the atmosphere temperature is \ + warmer than the ocean. This is not a situation in which the \ + thermodynamics are setup for right now." + end + new{FT, GT, DT}(grid, ocean, atmos, domain, floes) + end + + Model( + grid::GT, + ocean::Ocean{FT}, + atmos::Atmos{FT}, + domain::DT, + floes::StructArray{Floe{FT}}, + ) where { + FT<:AbstractFloat, + GT<:AbstractGrid{FT}, + DT<:Domain{ + FT, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + <:AbstractBoundary, + }, + } = + Model{FT, GT, DT}(grid, ocean, atmos, domain, floes) +end \ No newline at end of file diff --git a/src/simulation_components/oceans.jl b/src/simulation_components/oceans.jl new file mode 100644 index 0000000..0410029 --- /dev/null +++ b/src/simulation_components/oceans.jl @@ -0,0 +1,171 @@ +""" +Structs and functions used to define a Subzero oceans +""" + +""" + IceStressCell{FT<:AbstractFloat} + +Struct to collect stress from ice floes on ocean grid cells. One IceStressCell +corresponds to one grid cell. It holds a list of running totals of stress on the +cell, and a running list of the number of points making up those running totals. +Each element in the list corresponds to one floe, which is denoted in the +corresponding CellFloes matrix within the grid. +""" +struct IceStressCell{FT<:AbstractFloat} + τx::Vector{FT} + τy::Vector{FT} + npoints::Vector{Int} +end + +""" + IceStressCell{FT}() + +Constructs an IceStressCell object with empty lists for fields. +""" +IceStressCell{FT}() where {FT} = IceStressCell{FT}( + Vector{FT}(), + Vector{FT}(), + Vector{Int}() +) +""" + empty!(cell::IceStressCell) + +Empties the vectors within an IceStressCell +""" +function Base.empty!(cell::IceStressCell) + empty!(cell.τx) + empty!(cell.τy) + empty!(cell.npoints) + return +end + +""" + Ocean{FT<:AbstractFloat} + +Simulation ocean holding ocean values on the grid-scale with matricies of the +same size as the model's grid. The struct has the following fields: +- u is the ocean velocities in the x-direction for each grid cell +- v is the ocean velocities in the y-direction for each grid cell +- temp is the ocean temperature for each grid cell +- hflx_factor is a factor to calculate the ocean-atmosphere heat flux for a + cell in that grid cell by multiplying by its height +- si_frac is the fraction of area in each grid cell that is covered in sea-ice + +Ocean fields must all be matricies with dimensions equal to the number of grid +lines in the model's grid. +Note: If a periodic boundary is used in the domain, the last grid cell in that +direction will not be used as it is equivalent to the first grid cell. Thus, for +all of these fields, the first and last value in the x and/or y direction should +be equal if the east-west or north-south boundary pair are periodic +respectively. +""" +struct Ocean{FT<:AbstractFloat} + u::Matrix{FT} + v::Matrix{FT} + temp::Matrix{FT} + hflx_factor::Matrix{FT} + scells::Matrix{IceStressCell{FT}} + τx::Matrix{FT} + τy::Matrix{FT} + si_frac::Matrix{FT} + dissolved::Matrix{FT} + + function Ocean{FT}( + u, + v, + temp, + hflx, + scells, + τx, + τy, + si_frac, + dissolved, + ) where {FT <: AbstractFloat} + if !all(-3 .<= temp .<= 0) + @warn "Ocean temperatures are above the range for freezing. The \ + thermodynamics aren't currently setup for these conditions." + end + if !(size(u) == size(v) == size(τx) == size(τy)) + throw(ArgumentError("One or more of the ocean vector fields aren't \ + the same dimension.")) + end + if !(size(temp) == size(hflx) == size(si_frac) == size(dissolved)) + throw(ArgumentError("One or more of the ocean tracer fields aren't \ + the same dimension.")) + end + new{FT}(u, v, temp, hflx, scells, τx, τy, si_frac, dissolved) + end +end + +""" + Ocean(::Type{FT}, args...) + +A float type FT can be provided as the first argument of any Ocean constructor. +An Ocean of type FT will be created by passing all other arguments to the +correct constructor. +""" +Ocean(::Type{FT}, args...) where {FT <: AbstractFloat} = + Ocean{FT}(args...) + +""" + Ocean(args...) + +If a type isn't specified, Ocean will be of type Float64 and the correct +constructor will be called with all other arguments. +""" +Ocean(args...) = Ocean{Float64}(args...) + +""" + Ocean{FT}(u, v, temp) + +Construct model ocean. +Inputs: + u ocean x-velocity matrix with u for each grid line + v ocean y-velocity matrix with u for each grid line + temp temperature matrix with ocean/ice interface temperature for + each grid line +Output: + Model ocean with given velocity and temperature fields on each grid line. +""" +function Ocean{FT}( + u, + v, + temp, +) where {FT <: AbstractFloat} + Nx, Ny = size(temp) .- 1 + return Ocean{FT}( + u, + v, + temp, + zeros(FT, Nx + 1, Ny + 1), # heat flux + [IceStressCell{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], + zeros(FT, Nx + 1, Ny + 1), # x-stress + zeros(FT, Nx + 1, Ny + 1), # y-stress + zeros(FT, Nx + 1, Ny + 1), # sea ice fraction + zeros(FT, Nx + 1, Ny + 1), # dissolved + ) +end + +""" + Ocean{FT}(grid, u, v, temp) + +Construct model's ocean. +Inputs: + grid model grid + u ocean x-velocity for each grid line + v ocean y-velocity for each grid line + temp temperature at ocean/ice interface per grid cell +Output: + Ocean with constant velocity and temperature on each grid line. +""" +Ocean{FT}( + grid::AbstractGrid, + u, + v, + temp, +) where {FT <: AbstractFloat} = + Ocean{FT}( + fill(FT(u), grid.Nx + 1, grid.Ny + 1), + fill(FT(v), grid.Nx + 1, grid.Ny + 1), + fill(FT(temp), grid.Nx + 1, grid.Ny + 1), + ) diff --git a/src/simulation.jl b/src/simulation_components/simulation.jl similarity index 100% rename from src/simulation.jl rename to src/simulation_components/simulation.jl From 7603ce09535ec39c5ba7c161fc473675f37acd8a Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Mon, 3 Jul 2023 09:49:29 -0700 Subject: [PATCH 2/2] Recombine some files --- src/Subzero.jl | 8 +- src/simulation_components/atmos.jl | 66 ----- .../{domains.jl => domain_and_grid.jl} | 234 ++++++++++++++++-- src/simulation_components/grids.jl | 206 --------------- .../{oceans.jl => ocean_and_atmos.jl} | 65 ++++- 5 files changed, 288 insertions(+), 291 deletions(-) delete mode 100644 src/simulation_components/atmos.jl rename src/simulation_components/{domains.jl => domain_and_grid.jl} (73%) delete mode 100644 src/simulation_components/grids.jl rename src/simulation_components/{oceans.jl => ocean_and_atmos.jl} (72%) diff --git a/src/Subzero.jl b/src/Subzero.jl index e7aa915..d049cf3 100644 --- a/src/Subzero.jl +++ b/src/Subzero.jl @@ -92,9 +92,11 @@ const RingVec{T} = R where { } # Model -include("floe.jl") +include("simulation_components/floe.jl") include("floe_utils.jl") -include("model.jl") +include("simulation_components/domain_and_grid.jl") +include("simulation_components/ocean_and_atmos.jl") +include("simulation_components/model.jl") # Physical Processes include("physical_processes/update_floe.jl") include("physical_processes/coupling.jl") @@ -109,5 +111,5 @@ include("tools/compare_files.jl") include("logger.jl") # Simulation include("output.jl") -include("simulation.jl") +include("simulation_components/simulation.jl") end \ No newline at end of file diff --git a/src/simulation_components/atmos.jl b/src/simulation_components/atmos.jl deleted file mode 100644 index b1f0d7b..0000000 --- a/src/simulation_components/atmos.jl +++ /dev/null @@ -1,66 +0,0 @@ -""" -Structs and functions used to define a Subzero atmosphere -""" - - -""" -Atmos velocities in the x-direction (u) and y-direction (v). u and v should -match the size of the corresponding model grid so that there is one x and y -velocity value for each grid cell. Atmos also needs temperature at the -atmosphere/ice interface in each grid cell. Model cannot be constructed if size -of atmos fields and grid do not match. -""" -struct Atmos{FT<:AbstractFloat} - u::Matrix{FT} - v::Matrix{FT} - temp::Matrix{FT} - function Atmos{FT}(u, v, temp) where {FT <: AbstractFloat} - if !(size(u) == size(v)) - throw(ArgumentError("One or more of the atmosphere vector fields \ - aren't the same dimension.")) - end - new{FT}(u, v, temp) - end -end - -""" - Atmos(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any Atmos constructor. -An Atmos of type FT will be created by passing all other arguments to the -correct constructor. -""" -Atmos(::Type{FT}, args...) where {FT <: AbstractFloat} = - Atmos{FT}(args...) - -""" - Atmos(args...) - -If a type isn't specified, Atmos will be of type Float64 and the correct -constructor will be called with all other arguments. -""" -Atmos(args...) = Atmos{Float64}(args...) - -""" - Atmos{FT}(grid, u, v) - -Construct model atmosphere of type FT. -Inputs: - grid model's grid - u Atmos x-velocity for each grid cell - v Atmos y-velocity for each grid cell - temp temperature at atmopshere/ice interface per grid cell -Output: - Atmosphere of type FT with constant velocity and temperature over domain. -""" -Atmos{FT}( - grid::AbstractGrid, - u, - v, - temp, -) where {FT <: AbstractFloat} = - return Atmos{FT}( - fill(FT(u), grid.Nx + 1, grid.Ny + 1), - fill(FT(v), grid.Nx + 1, grid.Ny + 1), - fill(FT(temp), grid.Nx + 1, grid.Ny + 1), - ) diff --git a/src/simulation_components/domains.jl b/src/simulation_components/domain_and_grid.jl similarity index 73% rename from src/simulation_components/domains.jl rename to src/simulation_components/domain_and_grid.jl index 0ca64d9..993102b 100644 --- a/src/simulation_components/domains.jl +++ b/src/simulation_components/domain_and_grid.jl @@ -1,5 +1,5 @@ """ -Structs and functions used to define a Subzero domains +Structs and functions used to define a Subzero domain and grid """ """ @@ -43,7 +43,7 @@ rectangular domain. struct West<:AbstractDirection end """ - boundary_coords(grid::AbstractGrid, ::Type{North}) + boundary_coords(grid, ::Type{North}) Determine coordinates of northen-most boundary of domain if around the edge of the grid. @@ -59,7 +59,7 @@ Output: possibly too cautious. If boundary_coords methods are used for each direction, corners will be shared between adjacent boundaries. """ -function boundary_coords(grid::AbstractGrid, ::Type{North}) +function boundary_coords(grid, ::Type{North}) Δx = (grid.xf - grid.x0)/2 # Half of the grid in x Δy = (grid.yf - grid.y0)/2 # Half of the grid in y return grid.yf, # val @@ -71,7 +71,7 @@ function boundary_coords(grid::AbstractGrid, ::Type{North}) end """ - boundary_coords(grid::AbstractGrid, ::Type{South}) + boundary_coords(grid, ::Type{South}) Determine coordinates of southern-most boundary of domain if around the edge of the grid. @@ -82,7 +82,7 @@ Output: PolyVec of boundary coordinates. See documentation of North method of this function for more details. """ -function boundary_coords(grid::AbstractGrid, ::Type{South}) +function boundary_coords(grid, ::Type{South}) Δx = (grid.xf - grid.x0)/2 # Half of the grid in x Δy = (grid.yf - grid.y0)/2 # Half of the grid in y return grid.y0, # val @@ -94,7 +94,7 @@ function boundary_coords(grid::AbstractGrid, ::Type{South}) end """ - boundary_coords(grid::AbstractGrid, ::Type{East}) + boundary_coords(grid, ::Type{East}) Determine coordinates of eastern-most boundary of domain if around the edge of the grid. @@ -105,7 +105,7 @@ Output: PolyVec of boundary coordinates. See documentation of North method of this function for more details. """ -function boundary_coords(grid::AbstractGrid, ::Type{East}) +function boundary_coords(grid, ::Type{East}) Δx = (grid.xf - grid.x0)/2 # Half of the grid in x Δy = (grid.yf - grid.y0)/2 # Half of the grid in y return grid.xf, # val @@ -117,7 +117,7 @@ function boundary_coords(grid::AbstractGrid, ::Type{East}) end """ - boundary_coords(grid::AbstractGrid, ::Type{West}) + boundary_coords(grid, ::Type{West}) Determine coordinates of western-most boundary of domain if around the edge of the grid. @@ -128,7 +128,7 @@ Output: PolyVec of boundary coordinates. See documentation of North method of this function for more details. """ -function boundary_coords(grid::AbstractGrid, ::Type{West}) +function boundary_coords(grid, ::Type{West}) Δx = (grid.xf - grid.x0)/2 # Half of the grid in x Δy = (grid.yf - grid.y0)/2 # Half of the grid in y return grid.x0, # val @@ -208,7 +208,7 @@ OpenBoundary(::Type{D}, args...) where {D <: AbstractDirection} = OpenBoundary{D, Float64}(args...) """ - OpenBoundary{D, FT}(grid::AbstractGrid) + OpenBoundary{D, FT}(grid) Creates open boundary on the edge of the grid, and with the direction as a type. Edge is determined by direction. @@ -218,7 +218,7 @@ Outputs: Open Boundary on edge of grid given by direction and type. """ function OpenBoundary{D, FT}( - grid::AbstractGrid, + grid, ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) OpenBoundary{D, FT}( @@ -262,7 +262,7 @@ PeriodicBoundary(::Type{D}, args...) where {D <: AbstractDirection} = PeriodicBoundary{D, Float64}(args...) """ - PeriodicBoundary{D, FT}(grid::AbstractGrid) + PeriodicBoundary{D, FT}(grid) Creates periodic boundary on the edge of the grid, and with the direction as a type. Edge is determined by direction. @@ -272,7 +272,7 @@ Outputs: Periodic Boundary on edge of grid given by direction and type. """ function PeriodicBoundary{D, FT}( - grid::AbstractGrid, + grid, ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) PeriodicBoundary{D, FT}( @@ -328,7 +328,7 @@ Outputs: Collision Boundary on edge of grid given by direction. """ function CollisionBoundary{D, FT}( - grid::AbstractGrid, + grid, ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) CollisionBoundary{D, FT}( @@ -390,7 +390,7 @@ Note: should be negative. Else the velocity should be positive. """ function CompressionBoundary{D, FT}( - grid::AbstractGrid, + grid, velocity, ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) @@ -666,3 +666,207 @@ Domain( west, StructArray{TopographyElement{FT}}(undef, 0), ) + + """ + CellFloes{FT<:AbstractFloat} + + Struct that tracks which floes are within given cell, as well as the translation + vector needed to move floe each from current postion into cell if it is in cell + due to periodic boundaries. Each index in floeidx is the index of a floe within + the cell and the Δx and Δy with the same index are that floe's translation + vector. + Note: the floeidx is the index of grid cells centered on grid lines, not on the + grid cells defined by the regular, rectilinear grid. + """ + struct CellFloes{FT<:AbstractFloat} + floeidx::Vector{Int} + Δx::Vector{FT} + Δy::Vector{FT} + end + + """ + CellFloes{FT}() + + Constructs an CellFloes object with empty lists for fields. + """ + CellFloes{FT}() where {FT} = CellFloes{FT}( + Vector{Int}(), + Vector{FT}(), + Vector{FT}() + ) + """ + empty!(cell::CellFloes) + + Empties the vectors within a CellFloes object + """ + function Base.empty!(cell::CellFloes) + empty!(cell.floeidx) + empty!(cell.Δx) + empty!(cell.Δy) + return + end + + """ + AbstractGrid + + An abstract type for the grid that model will be simulated on. + Affects calculation on the grid. + """ + abstract type AbstractGrid{FT<:AbstractFloat} end + + """ + RegRectilinearGrid{FT<:AbstractFloat}<:AbstractGrid + + Tessellation of 2-dimensional Euclidean space into n-by-m congruent rectangles. + - Nx: number of grid cells in the x-direction + - Ny: number of grid cells in the y-direction + - x0: value of first x grid line + - xf: value of final x grid line + - y0: value of first y grid line + - yf: value of final y grid line + - Δx: grid cell width + - Δy: grid cell height + - floe_locations: matrix of CellFloes, which keep track of which floes are in + each cell + """ + struct RegRectilinearGrid{FT}<:AbstractGrid{FT} + Nx::Int + Ny::Int + x0::FT + xf::FT + y0::FT + yf::FT + Δx::FT + Δy::FT + floe_locations::Matrix{CellFloes{FT}} + + function RegRectilinearGrid{FT}( + Nx, + Ny, + x0, + xf, + y0, + yf, + Δx, + Δy, + floe_locations, + ) where {FT} + if size(floe_locations) != (Nx + 1, Ny + 1) + throw(ArgumentError("Floe location matrix needs to be the same \ + dimensions as grid lines.")) + end + if (xf - x0) / Δx != Nx + throw(ArgumentError("X-grid extent and grid cell width don't match \ + with number of grid cells in x-direction.")) + end + if (yf - y0) / Δy != Ny + throw(ArgumentError("Y-grid extent and grid cell height don't match \ + with number of grid cells in y-direction.")) + end + return new{FT}(Nx, Ny, x0, xf, y0, yf, Δx, Δy, floe_locations) + end + end + + """ + RegRectilinearGrid(::Type{FT}, args...) + + A float type FT can be provided as the first argument of any RegRectilinearGrid + constructor. A RegRectilinearGrid of type FT will be created by passing all + other arguments to the correct constructor. + """ + RegRectilinearGrid(::Type{FT}, args...) where {FT <: AbstractFloat} = + RegRectilinearGrid{FT}(args...) + + """ + RegRectilinearGrid(args...) + + If a type isn't specified, RegRectilinearGrid will be of type Float64 and the + correct constructor will be called with all other arguments. + """ + RegRectilinearGrid(args...) = RegRectilinearGrid{Float64}(args...) + + """ + RegRectilinearGrid{FT}( + xbounds::Tuple, + ybounds::Tuple, + Δx, + Δy, + ) + + Construct a RegRectilinearGrid for model given bounds for grid x and y and grid + cell dimensions in meters. + Inputs: + xbounds bound of grid x-direction in form (left, right) + ybounds bound of grid y-direction in form (bottom, top) + Δx length/height of grid cells in x-direction + Δy length/height of grid cells in y-direction + Output: + RegRectilinearGrid from lx to ux and height from ly to uy with grid squares + of size Δx by Δy + Warning: + If Δx doesn't evenly divide x length (lu-lx) or Δy doesn't evenly divide y + length (uy-ly) you won't get full size grid. The grid will be "trimmed" to + the nearest full grid square in both directions. + """ + function RegRectilinearGrid{FT}( + xbounds::Tuple, + ybounds::Tuple, + Δx, + Δy, + ) where {FT <: AbstractFloat} + Nx = floor(Int, (xbounds[2] - xbounds[1]) / Δx) + Ny = floor(Int, (ybounds[2] - ybounds[1]) / Δy) + return RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds[1], + xbounds[1] + Nx * Δx, + ybounds[1], + ybounds[1] + Ny * Δy, + Δx, + Δy, + [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], + ) + end + + get_grid_lines(g0, gf, Δg) = g0:Δg:gf + get_grid_centers(g0, gf, Δg) = (g0 + Δg/2):Δg:(gf - Δg/2) + + """ + RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds::Tuple{Real, Real}, + ybounds::Tuple{Real, Real}, + ) where {FT <: AbstractFloat} + + Construct a RegRectilinearGrid for model given bounds for grid x and y and the + number of grid cells in both the x and y direction. + Inputs: + Nx number of grid cells in the x-direction + Ny number of grid cells in the y-direction + xbounds bound of grid x-direction in form (left, right) + ybounds bound of grid y-direction in form (bottom, top) + Output: + RegRectilinearGrid with width and height determined by xbound and ybounds + and the number of grid cells in the x-direction and y-direction determined + by dims. + """ + RegRectilinearGrid{FT}( + Nx::Int, + Ny::Int, + xbounds::Tuple, + ybounds::Tuple, + ) where {FT <: AbstractFloat} = + RegRectilinearGrid{FT}( + Nx, + Ny, + xbounds[1], + xbounds[2], + ybounds[1], + ybounds[2], + (xbounds[2] - xbounds[1]) / Nx, + (ybounds[2] - ybounds[1]) / Ny, + [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], + ) + \ No newline at end of file diff --git a/src/simulation_components/grids.jl b/src/simulation_components/grids.jl deleted file mode 100644 index 2a60ab7..0000000 --- a/src/simulation_components/grids.jl +++ /dev/null @@ -1,206 +0,0 @@ -""" -Structs and functions used to define a Subzero grids -""" - -""" - CellFloes{FT<:AbstractFloat} - -Struct that tracks which floes are within given cell, as well as the translation -vector needed to move floe each from current postion into cell if it is in cell -due to periodic boundaries. Each index in floeidx is the index of a floe within -the cell and the Δx and Δy with the same index are that floe's translation -vector. -Note: the floeidx is the index of grid cells centered on grid lines, not on the -grid cells defined by the regular, rectilinear grid. -""" -struct CellFloes{FT<:AbstractFloat} - floeidx::Vector{Int} - Δx::Vector{FT} - Δy::Vector{FT} -end - -""" - CellFloes{FT}() - -Constructs an CellFloes object with empty lists for fields. -""" -CellFloes{FT}() where {FT} = CellFloes{FT}( - Vector{Int}(), - Vector{FT}(), - Vector{FT}() -) -""" - empty!(cell::CellFloes) - -Empties the vectors within a CellFloes object -""" -function Base.empty!(cell::CellFloes) - empty!(cell.floeidx) - empty!(cell.Δx) - empty!(cell.Δy) - return -end - -""" - AbstractGrid - -An abstract type for the grid that model will be simulated on. -Affects calculation on the grid. -""" -abstract type AbstractGrid{FT<:AbstractFloat} end - -""" - RegRectilinearGrid{FT<:AbstractFloat}<:AbstractGrid - -Tessellation of 2-dimensional Euclidean space into n-by-m congruent rectangles. -- Nx: number of grid cells in the x-direction -- Ny: number of grid cells in the y-direction -- x0: value of first x grid line -- xf: value of final x grid line -- y0: value of first y grid line -- yf: value of final y grid line -- Δx: grid cell width -- Δy: grid cell height -- floe_locations: matrix of CellFloes, which keep track of which floes are in - each cell -""" -struct RegRectilinearGrid{FT}<:AbstractGrid{FT} - Nx::Int - Ny::Int - x0::FT - xf::FT - y0::FT - yf::FT - Δx::FT - Δy::FT - floe_locations::Matrix{CellFloes{FT}} - - function RegRectilinearGrid{FT}( - Nx, - Ny, - x0, - xf, - y0, - yf, - Δx, - Δy, - floe_locations, - ) where {FT} - if size(floe_locations) != (Nx + 1, Ny + 1) - throw(ArgumentError("Floe location matrix needs to be the same \ - dimensions as grid lines.")) - end - if (xf - x0) / Δx != Nx - throw(ArgumentError("X-grid extent and grid cell width don't match \ - with number of grid cells in x-direction.")) - end - if (yf - y0) / Δy != Ny - throw(ArgumentError("Y-grid extent and grid cell height don't match \ - with number of grid cells in y-direction.")) - end - return new{FT}(Nx, Ny, x0, xf, y0, yf, Δx, Δy, floe_locations) - end -end - -""" - RegRectilinearGrid(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any RegRectilinearGrid -constructor. A RegRectilinearGrid of type FT will be created by passing all -other arguments to the correct constructor. -""" -RegRectilinearGrid(::Type{FT}, args...) where {FT <: AbstractFloat} = - RegRectilinearGrid{FT}(args...) - -""" - RegRectilinearGrid(args...) - -If a type isn't specified, RegRectilinearGrid will be of type Float64 and the -correct constructor will be called with all other arguments. -""" -RegRectilinearGrid(args...) = RegRectilinearGrid{Float64}(args...) - -""" - RegRectilinearGrid{FT}( - xbounds::Tuple, - ybounds::Tuple, - Δx, - Δy, - ) - -Construct a RegRectilinearGrid for model given bounds for grid x and y and grid -cell dimensions in meters. -Inputs: - xbounds bound of grid x-direction in form (left, right) - ybounds bound of grid y-direction in form (bottom, top) - Δx length/height of grid cells in x-direction - Δy length/height of grid cells in y-direction -Output: - RegRectilinearGrid from lx to ux and height from ly to uy with grid squares - of size Δx by Δy -Warning: - If Δx doesn't evenly divide x length (lu-lx) or Δy doesn't evenly divide y - length (uy-ly) you won't get full size grid. The grid will be "trimmed" to - the nearest full grid square in both directions. -""" -function RegRectilinearGrid{FT}( - xbounds::Tuple, - ybounds::Tuple, - Δx, - Δy, -) where {FT <: AbstractFloat} - Nx = floor(Int, (xbounds[2] - xbounds[1]) / Δx) - Ny = floor(Int, (ybounds[2] - ybounds[1]) / Δy) - return RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds[1], - xbounds[1] + Nx * Δx, - ybounds[1], - ybounds[1] + Ny * Δy, - Δx, - Δy, - [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], - ) -end - -get_grid_lines(g0, gf, Δg) = g0:Δg:gf -get_grid_centers(g0, gf, Δg) = (g0 + Δg/2):Δg:(gf - Δg/2) - -""" - RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds::Tuple{Real, Real}, - ybounds::Tuple{Real, Real}, - ) where {FT <: AbstractFloat} - -Construct a RegRectilinearGrid for model given bounds for grid x and y and the -number of grid cells in both the x and y direction. -Inputs: - Nx number of grid cells in the x-direction - Ny number of grid cells in the y-direction - xbounds bound of grid x-direction in form (left, right) - ybounds bound of grid y-direction in form (bottom, top) -Output: - RegRectilinearGrid with width and height determined by xbound and ybounds - and the number of grid cells in the x-direction and y-direction determined - by dims. -""" -RegRectilinearGrid{FT}( - Nx::Int, - Ny::Int, - xbounds::Tuple, - ybounds::Tuple, -) where {FT <: AbstractFloat} = - RegRectilinearGrid{FT}( - Nx, - Ny, - xbounds[1], - xbounds[2], - ybounds[1], - ybounds[2], - (xbounds[2] - xbounds[1]) / Nx, - (ybounds[2] - ybounds[1]) / Ny, - [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], - ) diff --git a/src/simulation_components/oceans.jl b/src/simulation_components/ocean_and_atmos.jl similarity index 72% rename from src/simulation_components/oceans.jl rename to src/simulation_components/ocean_and_atmos.jl index 0410029..267aa47 100644 --- a/src/simulation_components/oceans.jl +++ b/src/simulation_components/ocean_and_atmos.jl @@ -1,5 +1,5 @@ """ -Structs and functions used to define a Subzero oceans +Structs and functions used to define a Subzero ocean and atmosphere """ """ @@ -169,3 +169,66 @@ Ocean{FT}( fill(FT(v), grid.Nx + 1, grid.Ny + 1), fill(FT(temp), grid.Nx + 1, grid.Ny + 1), ) + + +""" +Atmos velocities in the x-direction (u) and y-direction (v). u and v should +match the size of the corresponding model grid so that there is one x and y +velocity value for each grid cell. Atmos also needs temperature at the +atmosphere/ice interface in each grid cell. Model cannot be constructed if size +of atmos fields and grid do not match. +""" +struct Atmos{FT<:AbstractFloat} + u::Matrix{FT} + v::Matrix{FT} + temp::Matrix{FT} + function Atmos{FT}(u, v, temp) where {FT <: AbstractFloat} + if !(size(u) == size(v)) + throw(ArgumentError("One or more of the atmosphere vector fields \ + aren't the same dimension.")) + end + new{FT}(u, v, temp) + end +end + +""" + Atmos(::Type{FT}, args...) + +A float type FT can be provided as the first argument of any Atmos constructor. +An Atmos of type FT will be created by passing all other arguments to the +correct constructor. +""" +Atmos(::Type{FT}, args...) where {FT <: AbstractFloat} = + Atmos{FT}(args...) + +""" + Atmos(args...) + +If a type isn't specified, Atmos will be of type Float64 and the correct +constructor will be called with all other arguments. +""" +Atmos(args...) = Atmos{Float64}(args...) + +""" + Atmos{FT}(grid, u, v) + +Construct model atmosphere of type FT. +Inputs: + grid model's grid + u Atmos x-velocity for each grid cell + v Atmos y-velocity for each grid cell + temp temperature at atmopshere/ice interface per grid cell +Output: + Atmosphere of type FT with constant velocity and temperature over domain. +""" +Atmos{FT}( + grid::AbstractGrid, + u, + v, + temp, +) where {FT <: AbstractFloat} = + return Atmos{FT}( + fill(FT(u), grid.Nx + 1, grid.Ny + 1), + fill(FT(v), grid.Nx + 1, grid.Ny + 1), + fill(FT(temp), grid.Nx + 1, grid.Ny + 1), + )