diff --git a/Project.toml b/Project.toml index c8159dc..3daffa0 100644 --- a/Project.toml +++ b/Project.toml @@ -50,6 +50,7 @@ julia = "1.9" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" [targets] -test = ["Test"] +test = ["LibGEOS", "Test"] diff --git a/README.md b/README.md index a35be10..e461e20 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ ## _Fast and Flexible Sea Ice Dynamics_ Subzero.jl is a native [Julia](https://julialang.org/) discrete-element model (DEM) for exploring fine-scale sea ice dynamics, -reimplementing MATLAB model [SubZero by Manucharyan and Montemuro](https://doi.org/10.1029/2022MS003247). +reimplementing and enhancing MATLAB model [SubZero by Manucharyan and Montemuro](https://doi.org/10.1029/2022MS003247). - 🚀 Runs over **35 times faster** that original MATLAB model for title simulation! - 🧩 Modular simulation model makes it easy to **customize simulations**! diff --git a/docs/Project.toml b/docs/Project.toml index b7895e5..89078b0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,8 @@ [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +GeoInterfaceMakie = "0edc0954-3250-4c18-859d-ec71c1660c08" +GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" +LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" - -[compat] -Documenter = "1.5" -Literate = "2.19" diff --git a/docs/literate/tutorial.jl b/docs/literate/tutorial.jl index c7b68a4..1059d57 100644 --- a/docs/literate/tutorial.jl +++ b/docs/literate/tutorial.jl @@ -7,7 +7,125 @@ # ## Tutorial - copy-pasteable version -# ## Core steps of an Subzero.jl simulation +# ## Core ideas behind Subzero.jl simulations -using Subzero # bring package into scope \ No newline at end of file +# The very first step of running a Subzero simulation is to bring the package into scope. +using Subzero # bring Subzero into scope +using CairoMakie, GeoInterfaceMakie # bring plotting packages into scope +import GeoInterface as GI + + +# ## Creating a Grid + +# Each Subzero model requires a grid object. The grid object defines the grid for the ocean +# and atmosphere. Ocean and atmosphere vectors (like `u` and `v` velocities) and tracers +# (like temperature) are recorded on these grid points and grid lines. The grid points are +# then used for interpolation for coupling between the ice, ocean, and atmosphere. + +# All Subzero grid objects are concrete types of the abstract type [`AbstractRectilinearGrid`](@ref). +# Right now, all grid objects must be Rectilinear. Currently the only implemented concrete type +# is a [`RegRectilinearGrid`](@ref). If you are interested in implementing another type of +# grid, see the [developer documentation]("devdocs.md"). + +# Here, we will go ahead and create an instance of `RegRectilinearGrid`. We need to specify +# the grid endpoints and either the number of grid cells in both directions, or the size of +# the grid cells. Here we will specity the number of grid cells in the x-direction, `Nx`, and +# in the y-direction, `Ny`. + +grid = RegRectilinearGrid(; x0 = -1e5, xf = 1e5, y0 = 0.0, yf = 1e5, Nx = 20, Ny = 10) + +# We plot a dashed box around the grid so that we can see the that the grid matches the extent given. +# We also place tick-marks at the desired grid cell lengths. Finally, set the plot's aspect ration to `2` +# as the x-extent is two-times larger than the y-extent. +fig = Figure(); +Axis(fig[1, 1]; # set up axis tick marks to match grid cells + xticks = range(grid.x0, grid.xf, 5), xminorticks = IntervalsBetween(5), + xminorticksvisible = true, xminorgridvisible = true, + yticks = range(grid.y0, grid.yf, 3), yminorticks = IntervalsBetween(5), + yminorticksvisible = true, yminorgridvisible = true, +) +lines!( # plot boundary of grid with a dashed line + [grid.x0, grid.x0, grid.xf, grid.xf, grid.x0], # point x-values + [grid.y0, grid.yf, grid.yf, grid.y0, grid.y0]; # point y-values + linestyle = :dash, linewidth = 3.0) +# Resize grid to layout +colsize!(fig.layout, 1, Aspect(1, 2)) +resize_to_layout!(fig) +fig # display the figure + +# ## Creating Boundaries + +# Next, each Subzero.jl model needs a `Domain`. A `Domain` defines the region of the grid that +# the ice floes are allowed in, what happens to them when they reach the boundaries of that +# region, and if there is any topography in the model, along with the ice, in that region. + +# Similarly to the `grid` above, the `Domain` will be rectilinear, defined by four boundaries, +# one for each of the cardinal direction. You will be able to pass each of the cardinal +# directions ([`North`](@ref), [`South`](@ref), [`East`](@ref), and [`West`](@ref)), defined as +# types by Subzero, to the boundary constructors. Each boundary can have different behavior, allowing +# the user to create a wide variety of domain behavior. Right now, four types of `AbstractBoundaries` +# are implemented: [`OpenBoundary`](@ref), [`PeriodicBoundary`](@ref), [`CollisionBoundary`](@ref), and +# [`MovingBoundary`](@ref). + +# In this example, we will use two `CollisionBoundary` walls and two `PeriodicBoundary` walls +# to create a channel that the ice can infinitly flow through, from the east back to the west. +# In the north and the south, the ice will collide with the boundaries, as if there was shoreline +# preventing it from leaving the channel. + +# We will use the grid we made above to define the boundaries so that they exactly border the grid. + +north_bound = CollisionBoundary(North; grid) +south_bound = CollisionBoundary(South; grid) +east_bound = PeriodicBoundary(East; grid) +west_bound = PeriodicBoundary(West; grid) + +# If we plot the polygons that are created within each boundary object, we can see that they border +# the grid. These polygons are how we well when the ice floes are interacting with each of the +# boundaries. We can also see that the boundaries overlap in the corners to ensure there is +# a solid border around the grid. The `PeriodicBoundary` elements are in purple while the +# `CollisionBoundary` elements are in teal. + +poly!( # plot each of the boundaries with a 50% transparent color so we can see the overlap + [north_bound.poly, south_bound.poly, east_bound.poly, west_bound.poly]; + color = [(:purple, 0.5), (:purple, 0.5), (:teal, 0.5), (:teal, 0.5)], +) +fig # display the figure + +# ## Creating Topography +# We then have the option to add in a [`TopographyField`](@ref Subzero.TopographyField), which is a collection of +# [`TopographyElement`](@ref)s. If we want to add in topography field, we can create one using +# the [`initialize_topography_field`](@ref) function. Here we will create two islands in the +# channel. For simplcity, both will be triangles. I create the polygons that define the shape +# of each island using [`GeoInterface`](https://github.com/JuliaGeo/GeoInterface.jl) and +# defining the points with tuples: +island1 = GI.Polygon([[(-6e4, 7.5e4), (-4e4, 5e4), (-2.5e4, 7e4), (-6e4, 7.5e4)]]) +island2 = GI.Polygon([[(5e4, 2.5e4), (5e4, 5.5e4), (7.5e4, 3e4), (5e4, 2.5e4)]]) +topo_list = [island1, island2] + +# We can then pass these to `initialize_topography_field` with the `polys` keyword. We could +# also have defined them just by their coordinates and passed in the coordiantes by the `coords` +# keyword. +topo_field = initialize_topography_field(; polys = topo_list) + +# We can now plot the topography within the domain. +topo_color = RGBf(115/255, 93/255, 55/255) # brown color for topography +poly!(topo_field.poly; color = topo_color) # plot the topography +fig # display the figure + + +# ## Creating a Domian +# We now have all of the pieces needed to create a [`Domain`](@ref). We will combine the four +# (4) boundaries we created, and the topography, into one `Domain` object. The collection of +# boundaries and topography define where the floes can and cannot go in the simulation and add +# in boundary behavior.We can do that as follows: + +domain = Domain(; north = north_bound, south = south_bound, east = east_bound, west = west_bound, topography = topo_field) + +# !!! note +# We could have skipped adding the topography to the domain. That is an optional +# keyword and an empty topography field will be automatically created if the user does not +# provide one. + +# At this point, we have already plotted all of the `Domain` objects, so we will move in to adding +# environmental forcing from the ocean and the atmosphere. diff --git a/docs/make.jl b/docs/make.jl index 8307c0c..73cfba3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,12 +1,11 @@ using Documenter, Literate - using Subzero +import GeoInterface as GI +import GeometryOps as GO +import LibGEOS as LG + # Converting any files in the literate folder to markdown -# Literate.markdown( -# joinpath(@__DIR__, "src", "tutorial.jl"), joinpath(@__DIR__, "src"); -# credit = false -# ) LITERATE_INPUT = joinpath(@__DIR__, "literate") LITERATE_OUTPUT = joinpath(@__DIR__, "src") @@ -21,6 +20,7 @@ for (root, _, files) ∈ walkdir(LITERATE_INPUT), file ∈ files Literate.markdown(ipath, opath) end +# Documentation formatting format = Documenter.HTML(; repolink = "https://github.com/Caltech-OCTO/Subzero.jl", canonical = "https://Caltech-OCTO.github.io/SubzeroDocumentation/stable", @@ -28,6 +28,9 @@ format = Documenter.HTML(; size_threshold = 5*10^5, # 500 KiB ) +# Metadata from doc tests +DocMeta.setdocmeta!(Subzero, :DocTestSetup, :(using Subzero); recursive=true) + makedocs(; modules=[Subzero], authors="Skylar Gering and contributers", diff --git a/docs/src/api.md b/docs/src/api.md index 8bf05be..5389053 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -7,12 +7,46 @@ CurrentModule = Subzero !!! warning This page is still very much WIP! -Documentation for [Subzero](https://github.com/Caltech-OCTO/Subzero.jl)'s full API (only for reference!). +## Grids -```@index +```@docs +AbstractRectilinearGrid +RegRectilinearGrid ``` -## All methods -```@autodocs -Modules = [Subzero] -``` \ No newline at end of file +## Directions +```@docs +AbstractDirection +North +South +East +West +``` +## Boundaries +```@docs +AbstractBoundary +OpenBoundary +PeriodicBoundary +CollisionBoundary +MovingBoundary +``` +## Topography +```@docs +TopographyElement +initialize_topography_field +``` +## Domain +```@docs +Domain +``` + +### Model +```@docs +Model +``` + +## Developer-Used Types +```@docs +CellFloes +TopographyField +``` diff --git a/docs/src/assets/title.jl b/docs/src/assets/title.jl index a11e6f4..7cfa56a 100644 --- a/docs/src/assets/title.jl +++ b/docs/src/assets/title.jl @@ -10,7 +10,7 @@ const Δt = 20 const nΔt = 15000 # Model instantiation -grid = RegRectilinearGrid((0.0, Lx), (0.0, Ly), Δgrid, Δgrid) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) uvels = repeat( [range(0, 0.2, length = 21); range(0.2, 0, length = 20)], @@ -24,11 +24,11 @@ ocean = Ocean( atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) # Floe creation floe_arr = initialize_floe_field(FT, 400, [0.85], domain, hmean, Δh; rng = Xoshiro(1)) diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 0ef3bce..918628f 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -10,12 +10,153 @@ simulatuion. ## Tutorial - copy-pasteable version -## Core steps of an Subzero.jl simulation +## Core ideas behind Subzero.jl simulations + +The very first step of running a Subzero simulation is to bring the package into scope. + +````@example tutorial +using Subzero # bring Subzero into scope +using CairoMakie, GeoInterfaceMakie # bring plotting packages into scope +import GeoInterface as GI +```` + +## Creating a Grid + +Each Subzero model requires a grid object. The grid object defines the grid for the ocean +and atmosphere. Ocean and atmosphere vectors (like `u` and `v` velocities) and tracers +(like temperature) are recorded on these grid points and grid lines. The grid points are +then used for interpolation for coupling between the ice, ocean, and atmosphere. + +All Subzero grid objects are concrete types of the abstract type [`AbstractRectilinearGrid`](@ref). +Right now, all grid objects must be Rectilinear. Currently the only implemented concrete type +is a [`RegRectilinearGrid`](@ref). If you are interested in implementing another type of +grid, see the [developer documentation]("devdocs.md"). + +Here, we will go ahead and create an instance of `RegRectilinearGrid`. We need to specify +the grid endpoints and either the number of grid cells in both directions, or the size of +the grid cells. Here we will specity the number of grid cells in the x-direction, `Nx`, and +in the y-direction, `Ny`. + +````@example tutorial +grid = RegRectilinearGrid(; x0 = -1e5, xf = 1e5, y0 = 0.0, yf = 1e5, Nx = 20, Ny = 10) +```` + +We plot a dashed box around the grid so that we can see the that the grid matches the extent given. +We also place tick-marks at the desired grid cell lengths. Finally, set the plot's aspect ration to `2` +as the x-extent is two-times larger than the y-extent. + +````@example tutorial +fig = Figure(); +Axis(fig[1, 1]; # set up axis tick marks to match grid cells + xticks = range(grid.x0, grid.xf, 5), xminorticks = IntervalsBetween(5), + xminorticksvisible = true, xminorgridvisible = true, + yticks = range(grid.y0, grid.yf, 3), yminorticks = IntervalsBetween(5), + yminorticksvisible = true, yminorgridvisible = true, +) +lines!( # plot boundary of grid with a dashed line + [grid.x0, grid.x0, grid.xf, grid.xf, grid.x0], # point x-values + [grid.y0, grid.yf, grid.yf, grid.y0, grid.y0]; # point y-values + linestyle = :dash, linewidth = 3.0) +```` + +Resize grid to layout + +````@example tutorial +colsize!(fig.layout, 1, Aspect(1, 2)) +resize_to_layout!(fig) +fig # display the figure +```` + +## Creating Boundaries + +Next, each Subzero.jl model needs a `Domain`. A `Domain` defines the region of the grid that +the ice floes are allowed in, what happens to them when they reach the boundaries of that +region, and if there is any topography in the model, along with the ice, in that region. + +Similarly to the `grid` above, the `Domain` will be rectilinear, defined by four boundaries, +one for each of the cardinal direction. You will be able to pass each of the cardinal +directions ([`North`](@ref), [`South`](@ref), [`East`](@ref), and [`West`](@ref)), defined as +types by Subzero, to the boundary constructors. Each boundary can have different behavior, allowing +the user to create a wide variety of domain behavior. Right now, four types of `AbstractBoundaries` +are implemented: [`OpenBoundary`](@ref), [`PeriodicBoundary`](@ref), [`CollisionBoundary`](@ref), and +[`MovingBoundary`](@ref). + +In this example, we will use two `CollisionBoundary` walls and two `PeriodicBoundary` walls +to create a channel that the ice can infinitly flow through, from the east back to the west. +In the north and the south, the ice will collide with the boundaries, as if there was shoreline +preventing it from leaving the channel. + +We will use the grid we made above to define the boundaries so that they exactly border the grid. + +````@example tutorial +north_bound = CollisionBoundary(North; grid) +south_bound = CollisionBoundary(South; grid) +east_bound = PeriodicBoundary(East; grid) +west_bound = PeriodicBoundary(West; grid) +```` + +If we plot the polygons that are created within each boundary object, we can see that they border +the grid. These polygons are how we well when the ice floes are interacting with each of the +boundaries. We can also see that the boundaries overlap in the corners to ensure there is +a solid border around the grid. The `PeriodicBoundary` elements are in purple while the +`CollisionBoundary` elements are in teal. + +````@example tutorial +poly!( # plot each of the boundaries with a 50% transparent color so we can see the overlap + [north_bound.poly, south_bound.poly, east_bound.poly, west_bound.poly]; + color = [(:purple, 0.5), (:purple, 0.5), (:teal, 0.5), (:teal, 0.5)], +) +fig # display the figure +```` + +## Creating Topography +We then have the option to add in a [`TopographyField`](@ref Subzero.TopographyField), which is a collection of +[`TopographyElement`](@ref)s. If we want to add in topography field, we can create one using +the [`initialize_topography_field`](@ref) function. Here we will create two islands in the +channel. For simplcity, both will be triangles. I create the polygons that define the shape +of each island using [`GeoInterface`](https://github.com/JuliaGeo/GeoInterface.jl) and +defining the points with tuples: ````@example tutorial -using Subzero # bring package into scope +island1 = GI.Polygon([[(-6e4, 7.5e4), (-4e4, 5e4), (-2.5e4, 7e4), (-6e4, 7.5e4)]]) +island2 = GI.Polygon([[(5e4, 2.5e4), (5e4, 5.5e4), (7.5e4, 3e4), (5e4, 2.5e4)]]) +topo_list = [island1, island2] ```` +We can then pass these to `initialize_topography_field` with the `polys` keyword. We could +also have defined them just by their coordinates and passed in the coordiantes by the `coords` +keyword. + +````@example tutorial +topo_field = initialize_topography_field(; polys = topo_list) +```` + +We can now plot the topography within the domain. + +````@example tutorial +topo_color = RGBf(115/255, 93/255, 55/255) # brown color for topography +poly!(topo_field.poly; color = topo_color) # plot the topography +fig # display the figure +```` + +## Creating a Domian +We now have all of the pieces needed to create a [`Domain`](@ref). We will combine the four +(4) boundaries we created, and the topography, into one `Domain` object. The collection of +boundaries and topography define where the floes can and cannot go in the simulation and add +in boundary behavior.We can do that as follows: + +````@example tutorial +domain = Domain(; north = north_bound, south = south_bound, east = east_bound, west = west_bound, topography = topo_field) +```` + +!!! note + We could have skipped adding the topography to the domain. That is an optional + keyword and an empty topography field will be automatically created if the user does not + provide one. + +At this point, we have already plotted all of the `Domain` objects, so we will move in to adding +environmental forcing from the ocean and the atmosphere. + --- *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* diff --git a/documentation.md b/documentation.md index 0e593d3..bebc64a 100644 --- a/documentation.md +++ b/documentation.md @@ -65,7 +65,7 @@ A grid object has the following fields: Nx, Ny, x0, xf, y0, yf, Δx, and Δy Ny is the number of rows and Ny is the number of columns within the grid as defined by each grid cell. x0 and xf are the minimum and maximum x-bounds, and y0 and yf are the minimum and maximum y-bounds. Finally, Δx and Δy are the number of grid cells in the x and y directions. -Additionally, RegRectilinearGrid is a concrete subtype of AbstractGrid. More concrete subtypes may be added in the future. +Additionally, RegRectilinearGrid is a concrete subtype of AbstractRectilinearGrid. More concrete subtypes may be added in the future. ### Ocean: The ocean here represents 2D vector fields of the surface layer of the ocean. It includes the following: u-velocity, v-velocity, temperature, and stresses in the x and y direction, the amount of area per cell covered in ice, and the `hflx_factor` per cell, which allows local heat flux calculations based on the difference in atmosphere and ocean temperature in each cell. If you want to run Subzero without coupling, the ocean will be a set of prescribed fields. If you couple Subzero to Oceananigans, Oceananigans will provide the velocity and temperature fields, and Subzero will provide oceananigans with stress fields from the ice and atmosphere on the top layer of the ocean. diff --git a/examples/converge_diverge_flow.jl b/examples/converge_diverge_flow.jl index 788f041..c34d50c 100644 --- a/examples/converge_diverge_flow.jl +++ b/examples/converge_diverge_flow.jl @@ -10,12 +10,8 @@ const Δh = 0.0 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0, Lx), - (0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) + uvels = repeat( hcat( range(0.1 ,0.6, step = 0.1)', @@ -27,12 +23,12 @@ ocean = Ocean(uvels, zeros(grid.Nx + 1, grid.Ny + 1), zeros(grid.Nx + 1, grid.Ny atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) # Floe creation floe_arr = initialize_floe_field( diff --git a/examples/example_sims.jl b/examples/example_sims.jl index 2af60d5..1dbdd2c 100644 --- a/examples/example_sims.jl +++ b/examples/example_sims.jl @@ -25,10 +25,10 @@ floe1 = Floe(floe1_poly, hmean, Δh, u = 1.0) floe_arr = StructArray([floe1]) # Floe hitting topography -topo = TopographyElement([[[0.5e4, 5e4], [0.5e4, 7e4], [1e4, 7e4], [1e4, 5e4], [0.5e4, 5e4]]]) +topo = TopographyElement(; coords = [[[0.5e4, 5e4], [0.5e4, 7e4], [1e4, 7e4], [1e4, 5e4], [0.5e4, 5e4]]]) topo_arr = StructVector([topo for i in 1:1]) -domain = Subzero.Domain(nboundary, sboundary, eboundary, wboundary, topo_arr) +domain = Subzero.Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary, topography = topo_arr) floe1_poly = LG.Polygon([[[-1.75e4, 5e4], [-1.75e4, 7e4], [-1.25e4, 7e4], [-1.25e4, 5e4], [-1.75e4, 5e4]]]) diff --git a/examples/forcing_contained_floes.jl b/examples/forcing_contained_floes.jl index 6be8f05..1d1d577 100644 --- a/examples/forcing_contained_floes.jl +++ b/examples/forcing_contained_floes.jl @@ -10,12 +10,8 @@ const Δh = 0.0 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0.0, Lx), - (0.0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) + # Set ocean u velocities ocean_uvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1)) for i in CartesianIndices(ocean_uvels) @@ -48,12 +44,12 @@ ocean = Ocean( atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = OpenBoundary(North, grid) -sboundary = OpenBoundary(South, grid) -eboundary = OpenBoundary(East, grid) -wboundary = OpenBoundary(West, grid) +nboundary = OpenBoundary(North; grid) +sboundary = OpenBoundary(South; grid) +eboundary = OpenBoundary(East; grid) +wboundary = OpenBoundary(West; grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) floe_settings = FloeSettings( subfloe_point_generator = SubGridPointsGenerator(grid, 2), diff --git a/examples/many_floes.jl b/examples/many_floes.jl index 303ac01..5ab70e3 100644 --- a/examples/many_floes.jl +++ b/examples/many_floes.jl @@ -15,25 +15,19 @@ const coarse_nx = 10 const coarse_ny = 10 # Model instantiation -grid = RegRectilinearGrid( - FT, - (-Lx, Lx), - (-Ly, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(FT; x0 = -Lx, xf = Lx, y0 = -Ly, yf = Ly, Δx = Δgrid, Δy = Δgrid) ocean = Ocean(FT, grid, -0.2, 0.0, -1.0) atmos = Atmos(FT, grid, 0.0, 0.0, -3.0) # Domain creation - boundaries and topography -nboundary = OpenBoundary(FT, North, grid) -sboundary = OpenBoundary(FT, South, grid) -eboundary = OpenBoundary(FT, East, grid) -wboundary = OpenBoundary(FT, West, grid) +nboundary = OpenBoundary(North, FT; grid) +sboundary = OpenBoundary(South, FT; grid) +eboundary = OpenBoundary(East, FT; grid) +wboundary = OpenBoundary(West, FT; grid) -domain = Subzero.Domain(nboundary, sboundary, eboundary, wboundary) +domain = Subzero.Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) # Floe instantiation nfloes = 100 diff --git a/examples/moving_bounds.jl b/examples/moving_bounds.jl index 918ef8d..a6a01fd 100644 --- a/examples/moving_bounds.jl +++ b/examples/moving_bounds.jl @@ -10,22 +10,17 @@ const Δh = 0.125 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0.0, Lx), - (0.0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) ocean = Ocean(grid, 0.0, 0.0, 0.0) atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = MovingBoundary(North, grid, 0.0, -0.1) -sboundary = MovingBoundary(South, grid, 0.0, 0.1) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) +nboundary = MovingBoundary(North; grid, u = 0.0, v = -0.1) +sboundary = MovingBoundary(South; grid, u = 0.0, v = 0.1) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) # Floe creation floe_arr = initialize_floe_field( diff --git a/examples/restart_sim.jl b/examples/restart_sim.jl index cb18aa4..68178fb 100644 --- a/examples/restart_sim.jl +++ b/examples/restart_sim.jl @@ -14,12 +14,7 @@ const uomax = 2 dirs = [joinpath("output/restart_sim", "run_" * string(i)) for i in 1:n_part_sim] # Build initial model / simulation -grid = RegRectilinearGrid( - (0.0, L), - (0.0, L), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = L, y0 = 0.0, yf = L, Δx = Δgrid, Δy = Δgrid) ngrid = Int(L/Δgrid) + 1 ygrid = range(0,L,ngrid) @@ -36,11 +31,11 @@ ocean = Ocean( atmos = Atmos(FT, grid, 0.0, 0.0, 0.0) -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary,west = wboundary) floe_settings = FloeSettings(subfloe_point_generator = SubGridPointsGenerator(grid, 2)) floe_arr = initialize_floe_field( diff --git a/examples/shear_flow.jl b/examples/shear_flow.jl index 30bc86c..87fe96d 100644 --- a/examples/shear_flow.jl +++ b/examples/shear_flow.jl @@ -10,12 +10,8 @@ const Δh = 0.0 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0.0, Lx), - (0.0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) + uvels = repeat( [range(0, 0.5, length = 26); range(0.5, 0, length = 25)], outer = (1, 51), @@ -28,12 +24,12 @@ ocean = Ocean( atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary) coupling_settings = CouplingSettings( two_way_coupling_on = true, diff --git a/examples/simple_strait.jl b/examples/simple_strait.jl index fb9654e..5ffb33e 100644 --- a/examples/simple_strait.jl +++ b/examples/simple_strait.jl @@ -10,31 +10,23 @@ const Δh = 0.0 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0.0, Lx), - (0.0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) ocean = Ocean(grid, 0.0, -0.3, 0.0) atmos = Atmos(grid, 0.0, 0.0, 0.0) # Domain creation -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = CollisionBoundary(East, grid) -wboundary = CollisionBoundary(West, grid) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = CollisionBoundary(East; grid) +wboundary = CollisionBoundary(West; grid) island = [[[6e4, 4e4], [6e4, 4.5e4], [6.5e4, 4.5e4], [6.5e4, 4e4], [6e4, 4e4]]] topo1 = [[[0, 0.0], [0, 1e5], [2e4, 1e5], [3e4, 5e4], [2e4, 0], [0.0, 0.0]]] topo2 = [[[8e4, 0], [7e4, 5e4], [8e4, 1e5], [1e5, 1e5], [1e5, 0], [8e4, 0]]] -topo_arr = initialize_topography_field( - FT, - [island, topo1, topo2] -) +topo_arr = initialize_topography_field(FT; coords = [island, topo1, topo2]) -domain = Domain(nboundary, sboundary, eboundary, wboundary, topo_arr) +domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary, topography = topo_arr) coupling_settings = CouplingSettings( two_way_coupling_on = true, diff --git a/examples/test_run.jl b/examples/test_run.jl index 3ea4355..32a1ebe 100644 --- a/examples/test_run.jl +++ b/examples/test_run.jl @@ -48,19 +48,19 @@ grid = RegRectilinearGrid( 1e4, ) topo_coords = [[[5e4, 5e4], [5e4, 7e4], [7e4, 7e4], [7e4, 5e4], [5e4, 5e4]]] -collision_domain = Subzero.Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), - initialize_topography_field([topo_coords]) +collision_domain = Subzero.Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = CollisionBoundary(East; grid), + west = CollisionBoundary(West; grid), + topography = initialize_topography_field(; coords = [topo_coords]) ) -periodic_domain = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), +periodic_domain = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), ) consts = Constants() @@ -135,11 +135,11 @@ zonal_ocn = Ocean(grid, 0.5, 0.0, 0.0) zero_atmos = Atmos(grid, 0.5, 0.0, 0.0) -domain = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), +domain = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), ) # Floe instantiation diff --git a/examples/uniform_flow.jl b/examples/uniform_flow.jl index 5354220..9e08134 100644 --- a/examples/uniform_flow.jl +++ b/examples/uniform_flow.jl @@ -10,22 +10,17 @@ const Δh = 0.0 const Δt = 20 # Model instantiation -grid = RegRectilinearGrid( - (0.0, Lx), - (0.0, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid) ocean = Ocean(grid, 0.1, 0.0, 0.0) atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = PeriodicBoundary(North, grid) -sboundary = PeriodicBoundary(South, grid) -eboundary = PeriodicBoundary(East, grid) -wboundary = PeriodicBoundary(West, grid) +nboundary = PeriodicBoundary(North; grid) +sboundary = PeriodicBoundary(South; grid) +eboundary = PeriodicBoundary(East; grid) +wboundary = PeriodicBoundary(West; grid) -domain = Domain(nboundary, sboundary, eboundary, wboundary) +domain = (; north = nboundary, south = sboundary, east = eboundary, west = wboundary) # Floe creation floe_arr = initialize_floe_field( diff --git a/src/Subzero.jl b/src/Subzero.jl index e4a4ece..29c6590 100644 --- a/src/Subzero.jl +++ b/src/Subzero.jl @@ -8,25 +8,9 @@ module Subzero end Subzero export - AbstractGrid, - RegRectilinearGrid, Ocean, Atmos, - OpenBC, - PeriodicBC, - CollisionBC, - CompressionBC, - AbstractBoundary, - MovingBoundary, - CollisionBoundary, - PeriodicBoundary, - OpenBoundary, - North, - South, - East, - West, Domain, - TopographyElement, Floe, Constants, Model, @@ -51,7 +35,6 @@ export torque, overlap, initialize_floe_field, - initialize_topography_field, NoFracture, HiblerYieldCurve, MohrsCone, @@ -66,11 +49,11 @@ export OutputWriters, check_energy_momentum_conservation_julia, IceStressCell, - CellFloes, MonteCarloPointsGenerator, SubGridPointsGenerator import Base.@kwdef # this is being exported as of version 1.9 +import Base.show import GeometryOps as GO import GeometryOps.GeoInterface as GI import GeometryOps.GeometryBasics as GB @@ -102,17 +85,34 @@ const RingVec{T} = R where { R <: AbstractArray{V}, } -const Polys{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing} where T <: AbstractFloat +const Polys{T, V} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing} where T -Base.convert(::Type{Polys{Float32}}, p::Polys{Float64}) = GO.tuples(p, Float32) -Base.convert(::Type{Polys{Float64}}, p::Polys{Float32}) = GO.tuples(p, Float64) +Base.convert(::Type{Polys{Float32}}, p::Polys{<:Real}) = GO.tuples(p, Float32) +Base.convert(::Type{Polys{Float64}}, p::Polys{<:Real}) = GO.tuples(p, Float64) + +const FT_DEF = "`FT::Type{<:AbstractFloat}`: Float type used to run the simulation, either \ +`Float64` (default) or `Float32`." +const POLY_DEF = "`poly::Polys{FT}`: Polygon used to represent the shape of a floe or topography" +const POLY_LIST_DEF = "`polys::Vector{<:Polygon}`: list of polygons meant to represent a field \ +of floes or topography elements. Polygons can be any polygon type that supports GeoInterface." +const COORDS_LIST_DEF = "`coords::Vector{<:PolyVec}`: list of polygon coordinates meant to \ +represent a field of floes or topography elements. PolyVec refers to a Vector{Vector{<:Points}} \ +where the points can be tuples, vectors, or static vectors and the innermost vector refers to each ring of the polygon." +const CENTROID_DEF = "`centroid::Vector{FT}`: Two-element vector meant to represent the (x, y) \ +point that is the centroid of either a floe or topography" +const RMAX_DEF = "`rmax::FT`: Float length representing the maximum radius of a floe or topography \ +from the centroid to any given vertex" # Types include("simulation_components/stress_calculators.jl") # Model +include("simulation_components/grids.jl") +include("simulation_components/domain_components/abstract_domains.jl") +include("simulation_components/domain_components/boundaries.jl") +include("simulation_components/domain_components/topography.jl") +include("simulation_components/domain_components/domains.jl") include("simulation_components/floe.jl") include("floe_utils.jl") -include("simulation_components/domain_and_grid.jl") include("simulation_components/ocean_and_atmos.jl") include("simulation_components/model.jl") # Physical Processes diff --git a/src/output.jl b/src/output.jl index b50bf61..024497a 100644 --- a/src/output.jl +++ b/src/output.jl @@ -299,7 +299,7 @@ end GridOutputWriter{FT}( outputs::Vector{Symbol}, Δtout, - grid::AbstractGrid, + grid::AbstractRectilinearGrid, dims; dir = ".", filename = "gridded_data.nc", @@ -328,7 +328,7 @@ Output: """ function GridOutputWriter{FT}( Δtout, - grid::AbstractGrid, + grid::AbstractRectilinearGrid, dims; outputs::Vector{Symbol} = collect(get_known_grid_outputs()), dir = ".", diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl index 49028c1..8d99ee2 100644 --- a/src/physical_processes/collisions.jl +++ b/src/physical_processes/collisions.jl @@ -188,7 +188,7 @@ function calc_elastic_forces( end """ - get_velocity( + _get_velocity( floe, x, y, @@ -203,7 +203,7 @@ Outputs: u u velocity at point (x, y) assuming it is on given floe v v velocity at point (x, y) assuming it is on given floe """ -function get_velocity( +function _get_velocity( floe::FloeType{FT}, x::FT, y::FT, @@ -213,29 +213,6 @@ function get_velocity( return u, v end -""" - get_velocity( - element::AbstractDomainElement{FT}, - x, - y, - ) - -Get velocity, which is 0m/s by default, of a point on topography element or -boundary. -""" -get_velocity( - element::AbstractDomainElement{FT}, - _, - _, -) where {FT} = - FT(0), FT(0) - -get_velocity( - element::MovingBoundary, - _, - _, -) = element.u, element.v - """ calc_friction_forces( ifloe, @@ -278,8 +255,8 @@ function calc_friction_forces( px = fpoints[i, 1] py = fpoints[i, 2] nnorm = sqrt(normal[i, 1]^2 + normal[i, 2]^2) - iu, iv = get_velocity(ifloe, px, py) - ju, jv = get_velocity(jfloe, px, py) + iu, iv = _get_velocity(ifloe, px, py) + ju, jv = _get_velocity(jfloe, px, py) udiff = iu - ju vdiff = iv - jv @@ -490,112 +467,6 @@ function floe_domain_element_interaction!( return end -""" - normal_direction_correct!( - forces, - fpoints, - boundary::AbstractBoundary{North, <:AbstractFloat}, - ) - -Zero-out forces that point in direction not perpendicular to North boundary wall. -Inputs: - force normal forces on each of the n regions - greater than a minimum area - fpoint point force is applied on each of the n - regions greater than a minimum area - boundary domain's northern - boundary -Outputs: - None. All forces in the x direction set to 0 if the point the force is - applied to is in the northern boundary. -""" -function normal_direction_correct!( - forces::Matrix{FT}, - fpoints, - boundary::AbstractBoundary{North, <:AbstractFloat}, -) where {FT} - forces[fpoints[:, 2] .>= boundary.val, 1] .= FT(0.0) - return -end - -""" - normal_direction_correct!( - forces, - fpoints, - boundary::AbstractBoundary{South, <:AbstractFloat}, - ) - -Zero-out forces that point in direction not perpendicular to South boundary wall. -See normal_direction_correct! on northern wall for more information -""" -function normal_direction_correct!( - forces::Matrix{FT}, - fpoints, - boundary::AbstractBoundary{South, <:AbstractFloat}, -) where {FT} - forces[fpoints[:, 2] .<= boundary.val, 1] .= FT(0.0) - return - end - -""" - normal_direction_correct!( - forces, - fpoints, - boundary::AbstractBoundary{East, <:AbstractFloat}, - ) - -Zero-out forces that point in direction not perpendicular to East boundary wall. -See normal_direction_correct! on northern wall for more information -""" -function normal_direction_correct!( - forces::Matrix{FT}, - fpoints, - boundary::AbstractBoundary{East, <:AbstractFloat}, -) where {FT} - forces[fpoints[:, 1] .>= boundary.val, 2] .= FT(0.0) - return -end - -""" - normal_direction_correct!( - forces, - fpoints, - boundary::AbstractBoundary{<:AbstractFloat, West}, - ) - -Zero-out forces that point in direction not perpendicular to West boundary wall. -See normal_direction_correct! on northern wall for more information -""" -function normal_direction_correct!( - forces::Matrix{FT}, - fpoints, - boundary::AbstractBoundary{West, <:AbstractFloat}, -) where {FT} - forces[fpoints[:, 1] .<= boundary.val, 2] .= FT(0.0) - return -end - -""" - normal_direction_correct!( - forces, - fpoints, - ::TopographyElement, - ) - -No forces should be zero-ed out in collidions with topography elements. -Inputs: - None used. -Outputs: - None. -""" -function normal_direction_correct!( - forces, - fpoints, - ::TopographyElement, -) - return -end - """ floe_domain_element_interaction!( floe, @@ -663,7 +534,7 @@ function floe_domain_element_interaction!( region_areas, force_factor, ) - normal_direction_correct!(normal_forces, fpoints, element) + _normal_direction_correct!(normal_forces, fpoints, element) # Calculate frictional forces at each force point np = size(fpoints, 1) if np > 0 @@ -685,60 +556,6 @@ function floe_domain_element_interaction!( return end -""" - update_boundary!(args...) - -No updates to boundaries that aren't compression boundaries. -""" -function update_boundary!( - boundary::Union{OpenBoundary, CollisionBoundary, PeriodicBoundary}, - Δt, -) - return -end -""" - update_boundary!(boundary, Δt) - -Move North/South compression boundaries by given velocity. Update coords and val -fields to reflect new position. -Inputs: - boundary - domain compression boundary - Δt number of seconds in a timestep -Outputs: - None. Move boundary North or South depending on velocity. -""" -function update_boundary!( - boundary::MovingBoundary{D, FT}, - Δt, -) where {D <: Union{North, South}, FT <: AbstractFloat} - Δd = boundary.v * Δt - boundary.val += Δd - translate!(boundary.coords, zero(FT), Δd) - _translate_poly(FT, boundary.poly, zero(FT), Δd) -end -""" - update_boundary!(boundary, Δt) - -Move East/West compression boundaries by given velocity. Update coords and val -fields to reflect new position. -Inputs: - boundary - domain compression boundary - Δt number of seconds in a timestep -Outputs: - None. Move boundary East/West depending on velocity. -""" -function update_boundary!( - boundary::MovingBoundary{D, FT}, - Δt, -) where {D <: Union{East, West}, FT <: AbstractFloat} - Δd = boundary.u * Δt - boundary.val += Δd - translate!(boundary.coords, Δd, zero(FT)) - _translate_poly(FT, boundary.poly, Δd, zero(FT)) -end - """ update_boundaries!(domain) @@ -746,10 +563,11 @@ Update each boundary in the domain. For now, this simply means moving compression boundaries by their velocities. """ function update_boundaries!(domain, Δt) - update_boundary!(domain.north, Δt) - update_boundary!(domain.south, Δt) - update_boundary!(domain.east, Δt) - update_boundary!(domain.west, Δt) + _update_boundary!(domain.north, Δt) + _update_boundary!(domain.south, Δt) + _update_boundary!(domain.east, Δt) + _update_boundary!(domain.west, Δt) + return end """ @@ -889,8 +707,7 @@ potential_interaction( centroid2, rmax1, rmax2, -) = ((centroid1[1] - centroid2[1])^2 + (centroid1[2] - centroid2[2])^2) < - (rmax1 + rmax2)^2 +) = ((centroid1[1] - centroid2[1])^2 + (centroid1[2] - centroid2[2])^2) < (rmax1 + rmax2)^2 """ timestep_collisions!( diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index 7e1ffbe..aca00e5 100644 --- a/src/physical_processes/coupling.jl +++ b/src/physical_processes/coupling.jl @@ -609,7 +609,7 @@ Calculates subfloe point's cartesian coordiantes, polar coordiantes, velocity and index within the grid. Inputs: floe floe - grid model's grid + grid model's grid domain model's domain mc_cart pre-allocated nx2 matrix for floe's monte carlo point's cartesian coordinates where the first @@ -820,7 +820,7 @@ Inputs: grid - nx2 matrix of indices where the first column is the grid column index and the second column is the grid row index for cells centered on grid lines - grid model grid + grid model grid domain model domain atmos model atmosphere ocean model ocean @@ -921,7 +921,7 @@ Inputs: xmax center cell maxumum x value ymin center cell minimum y value ymax center cell maximum y value - grid model's grid + grid model's grid type of north or south boundary - periodic pair type of east or west boundary - periodic pair Output: @@ -958,7 +958,7 @@ Inputs: xmax center cell maxumum x value ymin center cell minimum y value ymax center cell maximum y value - grid + grid type of either north or south boundary - not a periodic pair type of either east or west boundary - periodic @@ -1004,7 +1004,7 @@ Inputs: xmax center cell maxumum x value ymin center cell minimum y value ymax center cell maximum y value - grid + grid type of either north or south boundary - periodic pair type of either east or west boundary - not a @@ -1050,7 +1050,7 @@ Inputs: xmax center cell maxumum x value ymin center cell minimum y value ymax center cell maximum y value - grid + grid type of either north or south boundary - not a periodic pair type of either east or west boundary - not a @@ -1406,7 +1406,7 @@ Inputs: centered on grid lines τx_ocn x-stress caused by ocean on point τy_ocn y-stress caused by ocean on point - grid model's grid + grid model's grid domain model's domain cell_floes matrix of CellFloes, one for each grid cell @@ -1472,7 +1472,7 @@ stress on each grid cell per-floe in grid cell is also recorded for use in calc_two_way_coupling! Inputs: floes model's floe list - grid model's grid + grid model's grid atmos model's atmosphere ocean model's ocean domain model's domain @@ -1604,7 +1604,7 @@ Calculate effects of ice and atmosphere on the ocean and update ocean stress fields and sea ice fraction. Inputs: floes model's floes - grid model's grid + grid model's grid atmos model's atmosphere ocean model's ocean domain model's domain diff --git a/src/physical_processes/simplification.jl b/src/physical_processes/simplification.jl index 458711c..a8a6835 100644 --- a/src/physical_processes/simplification.jl +++ b/src/physical_processes/simplification.jl @@ -268,7 +268,7 @@ Remove floes marked for removal and dissolve floes smaller than minimum floe area if the dissolve setting is on. Inputs: floes model's floes - grid model's grid + grid model's grid domain model's domain dissolved ocean's dissolved field floe_settings simulation's settings for making floes diff --git a/src/simulation_components/domain_and_grid.jl b/src/simulation_components/domain_and_grid.jl deleted file mode 100644 index 4abba07..0000000 --- a/src/simulation_components/domain_and_grid.jl +++ /dev/null @@ -1,892 +0,0 @@ -""" -Structs and functions used to define a Subzero domain and grid -""" - -""" - AbstractDirection - -An abstract type for the boundary cardinal directions within model domain. -Boundary direction will control behavior of sea ice floes at edges of domain. -""" -abstract type AbstractDirection end - -""" - North<:AbstractDirection - -A simple direction type representing if a boundary is the northern boundary in a -rectangular domain. -""" -struct North<:AbstractDirection end - -""" - South<:AbstractDirection - -A simple direction type representing if a boundary is the southern boundary in a -rectangular domain. -""" -struct South<:AbstractDirection end - -""" - East<:AbstractDirection - -A simple direction type representing if a boundary is the eastern boundary in a -rectangular domain. -""" -struct East<:AbstractDirection end - -""" - West<:AbstractDirection - -A simple direction type representing if a boundary is the western boundary in a -rectangular domain. -""" -struct West<:AbstractDirection end - -""" - boundary_coords(grid, ::Type{North}) - -Determine coordinates of northen-most boundary of domain if around the edge of -the grid. -Inputs: - grid model grid - boundary direction type -Output: - PolyVec of boundary coordinates. These coordinates describe a rectangle that - has a length 2-times the length of the grid in the x-direction, centered on - the grid so that there is a buffer of half of the grid on either side. The - height is half of the grid in the y-direction. This buffer prevents pieces - of floes from passing outside the boundary before the next timestep - - possibly too cautious. If boundary_coords methods are used for each - direction, corners will be shared between adjacent boundaries. -""" -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 - [[[grid.x0 - Δx, grid.yf], # coords - [grid.x0 - Δx, grid.yf + Δy], - [grid.xf + Δx, grid.yf + Δy], - [grid.xf + Δx, grid.yf], - [grid.x0 - Δx, grid.yf]]] -end - -""" - boundary_coords(grid, ::Type{South}) - -Determine coordinates of southern-most boundary of domain if around the edge of -the grid. -Inputs: - grid model grid - boundary direction type -Output: - PolyVec of boundary coordinates. See documentation of North method of this - function for more details. -""" -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 - [[[grid.x0 - Δx, grid.y0 - Δy], # coords - [grid.x0 - Δx, grid.y0], - [grid.xf + Δx, grid.y0], - [grid.xf + Δx, grid.y0 - Δy], - [grid.x0 - Δx, grid.y0 - Δy]]] -end - -""" - boundary_coords(grid, ::Type{East}) - -Determine coordinates of eastern-most boundary of domain if around the edge of -the grid. -Inputs: - grid model grid - boundary direction type -Output: - PolyVec of boundary coordinates. See documentation of North method of this - function for more details. -""" -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 - [[[grid.xf, grid.y0 - Δy], # coords - [grid.xf, grid.yf + Δy], - [grid.xf + Δx, grid.yf + Δy], - [grid.xf + Δx, grid.y0 - Δy], - [grid.xf, grid.y0 - Δy]]] -end - -""" - boundary_coords(grid, ::Type{West}) - -Determine coordinates of western-most boundary of domain if around the edge of -the grid. -Inputs: - grid model grid - boundary direction -Output: - PolyVec of boundary coordinates. See documentation of North method of this - function for more details. -""" -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 - [[[grid.x0 - Δx, grid.y0 - Δy], # coords - [grid.x0 - Δx, grid.yf + Δy], - [grid.x0, grid.yf + Δy], - [grid.x0, grid.y0 - Δy], - [grid.x0 - Δx, grid.y0 - Δy]]] -end - -""" - AbstractDomainElement{FT<:AbstractFloat} - -An abstract type for all of the element that create the shape of the domain: -the 4 boundary walls that make up the rectangular domain and the topography -within the domain. -""" -abstract type AbstractDomainElement{FT<:AbstractFloat} end - - -""" - AbstractBoundary{D<:AbstractDirection, FT}<:AbstractDomainElement{FT} - -An abstract type for the types of boundaries at the edges of the model domain. -Boundary types will control behavior of sea ice floes at edges of domain. -The direction given by type D denotes which edge of a domain this boundary could -be and type FT is the simulation float type (e.g. Float64 or Float32). - -Each boundary type has the coordinates of the boudnary as a field. These should -be shapes that completely seal the domain, and should overlap on the corners as -seen in the example below: - ________________ -|__|____val___|__| <- North coordinates include corners -| | | | -| | | | <- East and west coordinates ALSO include corners -| | | | -Each bounday type also has a field called "val" that holds value that defines -the line y = val or x = val (depending on boundary direction), such that if the -floe crosses that line it would be partially within the boundary. -""" -abstract type AbstractBoundary{ - D<:AbstractDirection, - FT, -}<:AbstractDomainElement{FT} end - -""" - OpenBoundary <: AbstractBoundary - -A sub-type of AbstractBoundary that allows a floe to pass out of the domain edge -without any effects on the floe. -""" -struct OpenBoundary{D, FT}<:AbstractBoundary{D, FT} - poly::Polys{FT} - coords::PolyVec{FT} - val::FT -end - -""" - OpenBoundary(::Type{FT}, ::Type{D}, args...) - -A float type FT can be provided as the first argument of any OpenBoundary -constructor. The second argument D is the directional type of the boundary. -An OpenBoundary of type FT and D will be created by passing all other arguments -to the correct constructor. -""" -OpenBoundary(::Type{FT}, ::Type{D}, args...) where { - FT <: AbstractFloat, - D <: AbstractDirection, -} = OpenBoundary{D, FT}(args...) - -""" - OpenBoundary(::Type{D}, args...) - -If a float type isn't specified, OpenBoundary will be of type Float64 and the -correct constructor will be called with all other arguments. -""" -OpenBoundary(::Type{D}, args...) where {D <: AbstractDirection} = - OpenBoundary{D, Float64}(args...) - -""" - 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. -Inputs: - grid model grid -Outputs: - Open Boundary on edge of grid given by direction and type. -""" -function OpenBoundary{D, FT}( - grid, -) where {D <: AbstractDirection, FT <: AbstractFloat} - val, coords = boundary_coords(grid, D) - OpenBoundary{D, FT}( - make_polygon(convert(PolyVec{FT}, coords)), - coords, - val, - ) -end - -""" - PeriodicBoundary <: AbstractBoundary - -A sub-type of AbstractBoundary that moves a floe from one side of the domain to -the opposite side of the domain when it crosses the boundary, bringing the floe -back into the domain. -""" -struct PeriodicBoundary{D, FT}<:AbstractBoundary{D, FT} - poly::Polys{FT} - coords::PolyVec{FT} - val::FT -end - -""" - PeriodicBoundary(::Type{FT}, ::Type{D}, args...) - -A float type FT can be provided as the first argument of any PeriodicBoundary -constructor. The second argument D is the directional type of the boundary. -A PeriodicBoundary of type FT and D will be created by passing all other -arguments to the correct constructor. -""" -PeriodicBoundary(::Type{FT}, ::Type{D}, args...) where { - FT <: AbstractFloat, - D <: AbstractDirection, -} = PeriodicBoundary{D, FT}(args...) - -""" - PeriodicBoundary(::Type{D}, args...) - -If a float type isn't specified, PeriodicBoundary will be of type Float64 and -the correct constructor will be called with all other arguments. -""" -PeriodicBoundary(::Type{D}, args...) where {D <: AbstractDirection} = - PeriodicBoundary{D, Float64}(args...) - -""" - 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. -Inputs: - grid model grid -Outputs: - Periodic Boundary on edge of grid given by direction and type. -""" -function PeriodicBoundary{D, FT}( - grid, -) where {D <: AbstractDirection, FT <: AbstractFloat} - val, coords = boundary_coords(grid, D) - PeriodicBoundary{D, FT}( - make_polygon(convert(PolyVec{FT}, coords)), - coords, - val, - ) -end - -""" - CollisionBoundary <: AbstractBoundary - -A sub-type of AbstractBoundary that stops a floe from exiting the domain by -having the floe collide with the boundary. The boundary acts as an immovable, -unbreakable ice floe in the collision. -""" -struct CollisionBoundary{D, FT}<:AbstractBoundary{D, FT} - poly::Polys{FT} - coords::PolyVec{FT} - val::FT -end - -""" - CollisionBoundary(::Type{FT}, ::Type{D}, args...) - -A float type FT can be provided as the first argument of any CollisionBoundary -constructor. The second argument D is the directional type of the boundary. -A CollisionBoundary of type FT and D will be created by passing all other -arguments to the correct constructor. -""" -CollisionBoundary(::Type{FT}, ::Type{D}, args...) where { - FT <: AbstractFloat, - D <: AbstractDirection, -} = CollisionBoundary{D, FT}(args...) - -""" - CollisionBoundary(::Type{D}, args...) - -If a float type isn't specified, CollisionBoundary will be of type Float64 and -the correct constructor will be called with all other arguments. -""" -CollisionBoundary(::Type{D}, args...) where {D <: AbstractDirection} = - CollisionBoundary{D, Float64}(args...) - - -""" - CollisionBoundary{D, FT}(grid) - -Creates collision boundary on the edge of the grid, and with the direction as a -type. Edge is determined by direction. -Inputs: - grid model grid - direction direction of boundary wall -Outputs: - Collision Boundary on edge of grid given by direction. -""" -function CollisionBoundary{D, FT}( - grid, -) where {D <: AbstractDirection, FT <: AbstractFloat} - val, coords = boundary_coords(grid, D) - CollisionBoundary{D, FT}( - make_polygon(convert(PolyVec{FT}, coords)), - coords, - val, - ) -end - -""" - MovingBoundary <: AbstractBoundary - -A sub-type of AbstractBoundary that creates a floe along the boundary that moves -towards the center of the domain at the given velocity, compressing the ice -within the domain. This subtype is a mutable struct so that the coordinates and -val can be changed as the boundary moves. The u and v velocities are in [m/s]. - -Note that with a u-velocity, east and west walls move towards the center of the domain, -providing compressive stress, and with a v-velocity, the bounday creates a shear stress by -incorporating the velocity into friction calculations but doesn't actually move. This means -that the boundaries cannot move at an angle, distorting the shape of the domain regardless -the the combination of u and v velocities. The same, but opposite is true for the north -and south walls. -""" -mutable struct MovingBoundary{D, FT}<:AbstractBoundary{D, FT} - poly::Polys{FT} - coords::PolyVec{FT} - val::FT - u::FT - v::FT -end - -""" - MovingBoundary(::Type{FT}, ::Type{D}, args...) - -A float type FT can be provided as the first argument of any MovingBoundary -constructor. The second argument D is the directional type of the boundary. -A MovingBoundary of type FT and D will be created by passing all other -arguments to the correct constructor. -""" -MovingBoundary(::Type{FT}, ::Type{D}, args...) where { - FT <: AbstractFloat, - D <: AbstractDirection, -} = MovingBoundary{D, FT}(args...) - -""" - MovingBoundary(::Type{D}, args...) - -If a float type isn't specified, MovingBoundary will be of type Float64 and -the correct constructor will be called with all other arguments. -""" -MovingBoundary(::Type{D}, args...) where {D <: AbstractDirection} = - MovingBoundary{D, Float64}(args...) - - -""" - MovingBoundary{D, FT}(grid, velocity) - -Creates compression boundary on the edge of the grid, and with the direction as -a type. -Edge is determined by direction. -Inputs: - grid model grid - u u velocity of boundary - v v velocity of boundary -Outputs: - MovingBoundary on edge of grid given by direction. -""" -function MovingBoundary{D, FT}( - grid, - u, - v, -) where {D <: AbstractDirection, FT <: AbstractFloat} - val, coords = boundary_coords(grid, D) - MovingBoundary{D, FT}( - make_polygon(convert(PolyVec{FT}, coords)), - coords, - val, - u, - v, - ) -end - -""" - NonPeriodicBoundary - -Union of all non-peridic boundary types to use as shorthand for dispatch. -""" -const NonPeriodicBoundary = Union{ - OpenBoundary, - CollisionBoundary, - MovingBoundary, -} - -""" - TopographyElement{FT}<:AbstractDomainElement{FT} - -Singular topographic element with coordinates field storing where the element is -within the grid. These are used to create the desired topography within the -simulation and will be treated as islands or coastline within the model -in that they will not move or break due to floe interactions, but they will -affect floes. -""" -struct TopographyElement{FT}<:AbstractDomainElement{FT} - poly::Polys{FT} - coords::PolyVec{FT} - centroid::Vector{FT} - rmax::FT - - function TopographyElement{FT}( - poly, - coords, # TODO: Can remove when not used anymore in other parts of the code - centroid, - rmax, - ) where {FT <: AbstractFloat} - if rmax <= 0 - throw(ArgumentError("Topography element maximum radius must be \ - positive and non-zero.")) - end - poly = GO.ClosedRing()(poly) - new{FT}(poly, find_poly_coords(poly), centroid, rmax) - end -end - -""" - TopographyElement(::Type{FT}, args...) - -A float type FT can be provided as the first argument of any TopographyElement -constructor. A TopographyElement of type FT will be created by passing all other -arguments to the correct constructor. -""" -TopographyElement(::Type{FT}, args...) where {FT <: AbstractFloat} = - TopographyElement{FT}(args...) - -""" - TopographyElement(args...) - -If a type isn't specified, TopographyElement will be of type Float64 and the -correct constructor will be called with all other arguments. -""" -TopographyElement(args...) = TopographyElement{Float64}(args...) - -""" - TopographyElement{FT}(poly) - -Constructor for topographic element with Polygon -Inputs: - poly -Output: - Topographic element of abstract float type FT -""" -function TopographyElement{FT}(poly::Polys) where {FT <: AbstractFloat} - rmholes!(poly) - centroid = collect(GO.centroid(poly)) # TODO: Remove collect once type is changed - coords = find_poly_coords(poly) - rmax = calc_max_radius(poly, centroid, FT) - return TopographyElement{FT}( - poly, - coords, - centroid, - rmax, - ) -end - -""" - TopographyElement{FT}(coords) - -Constructor for topographic element with PolyVec coordinates -Inputs: - coords -Output: - Topographic element of abstract float type FT -""" -TopographyElement{FT}(coords::PolyVec) where {FT <: AbstractFloat} = - TopographyElement{FT}(make_polygon(convert(PolyVec{FT}, coords))) - -""" - initialize_topography_field(args...) - -If a type isn't specified, the list of TopographyElements will each be of type -Float64 and the correct constructor will be called with all other arguments. -""" -initialize_topography_field(args...) = - initialize_topography_field(Float64, args...) - -""" - initialize_topography_field( - ::Type{FT}, - coords, - ) - -Create a field of topography from a list of polygon coordiantes. -Inputs: - Type{FT} Type for grid's numberical fields - - determines simulation run type - coords list of polygon coords to make into floes -Outputs: - topo_arr list of topography elements - created from given polygon coordinates -""" -function initialize_topography_field( - ::Type{FT}, - coords, -) where {FT <: AbstractFloat} - topo_multipoly = GO.DiffIntersectingPolygons()(GI.MultiPolygon(coords)) - topo_arr = StructArray{TopographyElement{FT}}(undef, GI.npolygon(topo_multipoly)) - for (i, p) in enumerate(GI.getpolygon(topo_multipoly)) - topo_arr[i] = TopographyElement{FT}(p) - end - return topo_arr -end - -""" - periodic_compat(b1, b2) - -Checks if two boundaries are compatible as a periodic pair. This is true if they -are both periodic, or if neither are periodic. Otherwise, it is false. -""" -function periodic_compat(::PeriodicBoundary, ::PeriodicBoundary) - return true -end - -function periodic_compat(::PeriodicBoundary, _) - return false -end - -function periodic_compat(_, ::PeriodicBoundary) - return false -end - -function periodic_compat(_, _) - return true -end - -""" -Domain that holds 4 Boundary elements, forming a rectangle bounding the model -during the simulation, and a list of topography elements. - -In order to create a Domain, three conditions need to be met. First, if needs to -be periodically compatible. This means that pairs of opposite boundaries both -need to be periodic if one of them is periodic. Next, the value in the north -boundary must be greater than the south boundary and the value in the east -boundary must be greater than the west in order to form a valid rectangle. - -Note: The code depends on the boundaries forming a rectangle oriented along the -cartesian grid. Other shapes/orientations are not supported at this time. -""" -struct Domain{ - FT<:AbstractFloat, - NB<:AbstractBoundary{North, FT}, - SB<:AbstractBoundary{South, FT}, - EB<:AbstractBoundary{East, FT}, - WB<:AbstractBoundary{West, FT}, - TT<:StructArray{<:TopographyElement{FT}}, -} - north::NB - south::SB - east::EB - west::WB - topography::TT - - function Domain{FT, NB, SB, EB, WB, TT}( - north::NB, - south::SB, - east::EB, - west::WB, - topography::TT, - ) where { - FT<:AbstractFloat, - NB<:AbstractBoundary{North, FT}, - SB<:AbstractBoundary{South, FT}, - EB<:AbstractBoundary{East, FT}, - WB<:AbstractBoundary{West, FT}, - TT<:StructArray{<:TopographyElement{FT}}, - } - if !periodic_compat(north, south) - throw(ArgumentError("North and south boundary walls are not \ - periodically compatable as only one of them is periodic.")) - elseif !periodic_compat(east, west) - throw(ArgumentError("East and west boundary walls are not \ - periodically compatable as only one of them is periodic.")) - elseif north.val < south.val - throw(ArgumentError("North boundary value is less than south \ - boundary value.")) - elseif east.val < west.val - throw(ArgumentError("East boundary value is less than west \ - boundary value.")) - end - new{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) - end - - Domain( - north::NB, - south::SB, - east::EB, - west::WB, - topography::TT, - ) where { - FT<:AbstractFloat, - NB<:AbstractBoundary{North, FT}, - SB<:AbstractBoundary{South, FT}, - EB<:AbstractBoundary{East, FT}, - WB<:AbstractBoundary{West, FT}, - TT<:StructArray{<:TopographyElement{FT}}, - } = - Domain{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) -end - -function get_domain_element(domain, idx) - if idx == -1 - return domain.north - elseif idx == -2 - return domain.south - elseif idx == -3 - return domain.east - elseif idx == -4 - return domain.west - else - topo_idx = -(idx + 4) - return get_floe(domain.topography, topo_idx) - end -end - -""" - Domain(north, south, east, west) - -Creates domain with empty list of topography and given boundaries. -Inputs: - north north boundary - south south boundary - east east boundary - west west boundary -""" -Domain( - north::NB, - south::SB, - east::EB, - west::WB, -) where { - FT<:AbstractFloat, - NB<:AbstractBoundary{North, FT}, - SB<:AbstractBoundary{South, FT}, - EB<:AbstractBoundary{East, FT}, - WB<:AbstractBoundary{West, FT}, -} = - Domain( - north, - south, - east, - 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)], - ) - diff --git a/src/simulation_components/domain_components/abstract_domains.jl b/src/simulation_components/domain_components/abstract_domains.jl new file mode 100644 index 0000000..91e5a26 --- /dev/null +++ b/src/simulation_components/domain_components/abstract_domains.jl @@ -0,0 +1,151 @@ +export AbstractDomainElement # All elements that make up the domain +export AbstractDirection # # Cardinal directions for directional boundaries --> see boundaries.jl +export AbstractBoundary # Used to create concrete boundary types --> see boundaries.jl + +# Abstract type definitions and general fallback functions for all subtypes +""" + YourDomainElement{FT} <: AbstractDomainElement{FT} + +An abstract type for all of the element that create the shape and behavior of the [`Domain`](@ref): +the 4 boundary walls that make up the rectangular domain and the topography within the `Domain`. + +When building a [`Domain`](@ref), the user will create four (4) boundary elements that are +concrete types of [`AbstractBoundary`](@ref), which is itself a subtype of +`AbstractDomainElement`, and a list of [`TopographyElement`](@ref) objects, which are also a +concrete subtype of `AbstractDomainElement`. + +## _API_ +If the user wanted to add more domain elements, other than boundaries and topography, they +would create a new subtype of `AbstractDomainElement` and would need to write methods for +the following functions: +- `_get_velocity(element::AbstractDomainElement{FT}, x::FT, y::FT)` +- `_normal_direction_correct!(forces::Matrix{FT}, fpoints{Matrix{FT}}, element::AbstractDomainElement{FT})` + +The `_get_velocity` function gets the velocity of a domain `element` at the point (`x`, `y`). +Given that currently implemented domain elements cannot rotate, the point at which the +velocity is measured is irrelevant. Furthermore, all domain elements but +[`MovingBoundary`](@ref) are stationary and thus have zero-velocity. This function is called +in the `calc_friction_forces` function. + +The `_normal_direction_correct` function updates the `forces` provided to it and zero's out +any forces not in the normal direction, if applicable. Right now, this is used to ensure +that the elastic collision with boundaries only produces a normal foces on floes. For other +domain elements (topography), this function does nothing. + +_Notes_: +- If the user is wants to create a new boundary type, they must also implement the +function in the [`AbstractBoundary`](@ref) API. +- If the user wanted new elements in the [`Domain`](@ref), the user would need to change the +`Domain` struct. +""" +abstract type AbstractDomainElement{FT<:AbstractFloat} end + +""" + YourDirection <: AbstractDirection + +Each domain within a Subzero.jl [`Model`](@ref) must have four (4) boundaries (subtypes +of [`AbstractBoundary`](@ref)) where each of these boundaries is parametrically typed by the +direction of the boundary. The user will first choose one of the four cardinal directions, +the subtypes of `AbstractDirection`: +- [`North`](@ref) +- [`South`](@ref) +- [`East`](@ref) +- [`West`](@ref) + +This abstract type is not meant to be extended by the user, unless the user wants to move +away from a rectangular domain with assigned cardinal directions for each wall. This would +a more major redesign and the user should check out the [developer documentation]("devdocs.md"). + +## _API +- `_boundary_info_from_extent(D::Type{AbstractDirection}, FT::Type{AbstractFloat}, x0, xf, y0, yf)` + +The `_boundary_info_from_extent` function finds the values for the `poly` and `val` fields of +a subtype of [`AbstractBoundary`](@ref) of direction `D` given a region of the extent defined +by `x0`, `xf`, `y0`, and `yf` given that these are the minimum and maximum `x` and `y` values +of the region desired. The returned extent will have values of type `FT`. +""" +abstract type AbstractDirection end + +""" + AbstractBoundary{D, FT} <: AbstractDomainElement{FT} + +When running a Subzero simulation, each simulation has a rectangular domain that contains +the ice pack. This domain is made up for four boundaries, one in each of the cardinal directions: +[`North`](@ref), [`South`](@ref), [`East`](@ref), and [`West`](@ref). The user will pick +four boundary types, each concrete subtypes of `AbstractBoundary` and assign each of them a +direction, `D` and a float type `FT`, by typing them +[parametrically]("https://docs.julialang.org/en/v1/manual/types/#Parametric-Types"). + +The user can pick from the currently implemented concrete subtypes of `AbstractBoundary` +([`OpenBoundary`](@ref), [`CollisionBoundary`](@ref), [`PeriodicBoundary`](@ref), and +[`MovingBoundary`](@ref)) or implement their own. + +For now, each boundary wall is assumed to be a rectangle. For correct behavior, the corners +of each of the boundary rectangles should overlap to make sure that no floes reach outside +of the domain without touching a boundary. These rectangles should be represented as +polygons with the `poly` field within each boundary. Furthermore, each boundary should have +a `val` field that represents either the line `y = val` (for `North` and `South` walls) +or `x = val` (for `East` and `West` walls) that denotes the line that marks the "inner-most" +edge of the domain. If an ice floe passes over the boundary's `val` line, the floe interacts +with the boundary. + +```text + ________________ +|__|____val___|__| <- North coordinates include corners +| | | | +| | | | <- East and west coordinates ALSO include corners +| | | | +``` + +For ease of use, each boundary should have a constructor option where the user can provide +an [`AbstractRectilinearGrid`](@ref) concrete type to define a boundary where `val` is at +the edge of the grid so that the boundaries form a border right around the grid. + +## _API_ +In addition to the below function, any new `AbstractBoundary` subtype must also implement +[`AbstractDomainElement`](@ref) API functions. + +- `_update_boundary!(boundary::AbstractBoundary, Δt::Int)` +- `_periodic_compat(boundary1::AbstractBoundary, boundary2::AbstractBoundary)` + +The `_update_boundary!` function updates a boundary's position (by changing the `val` and +`poly` fields) at every timestep. Currently, only `MovingBoundary` elements are updated, and +their update depends on the length of the simulation's timestep, `Δt`. + +The `_periodic_compat` function determines if a pair of functions are "periodically compatible", +that is, could they form a pair on opposite edges of a domain to create a set of functioning +periodic boundaries. It is called when creating a [`Domain`](@ref). +""" +abstract type AbstractBoundary{ + D<:AbstractDirection, + FT, +}<:AbstractDomainElement{FT} end + +# Helper function for boundary pretty printing for boundaries with `poly` and `val` fields +function show_boundary_poly_val_strings(boundary::AbstractBoundary) + points = join(Set(GI.getpoint(boundary.poly))|>collect, ", ") + points_summary = "polygon points are defined by the following set: $points" + val_summary = "val is $(boundary.val)" + return points_summary, val_summary +end + +# Constants used in documentation for sub-type fields and keywords +const StaticQuadrilateral{FT} = GI.Polygon{false,false, SA.SVector{1, GI.LinearRing{false, false, SA.SVector{5, Tuple{FT, FT}}, Nothing, Nothing}},Nothing,Nothing} where FT + +const D_DEF = "`D::Type{<:AbstractDirection}`: subtype of AbstractDirection used to represent \\ +if a boundary is the North, South, East, or West boundary of a simulation" + +const BOUNDARY_POLY_DEF = "`poly::StaticQuadrilateral{FT}`: rectangular polygon representing \\ +the shape and location of the boundary on a given timestep with points of float type `FT`." + +const VAL_DEF = "`val::FT`: number of float type `FT` representing either the line `y = val` \\ +(for `North` and `South` walls) or `x = val` (for `East` and `West` walls) that denotes the \\ +line marking the inner-most edge of the domain." + +const X0_DEF = "`x0::FT`: minimum x-value of domain extent (region floes will be contained within)" + +const XF_DEF = "`xf::FT`: maximum x-value of domain extent (region floes will be contained within)" + +const Y0_DEF = "`y0::FT`: minimum y-value of domain extent (region floes will be contained within)" + +const YF_DEF = "`yf::FT`: maximum y-value of domain extent (region floes will be contained within)" diff --git a/src/simulation_components/domain_components/boundaries.jl b/src/simulation_components/domain_components/boundaries.jl new file mode 100644 index 0000000..18f71f9 --- /dev/null +++ b/src/simulation_components/domain_components/boundaries.jl @@ -0,0 +1,561 @@ +# Cardinal directions +export North, South, East, West +# Boundary of the model - types for each of the 4 walls (north, south, east, west) +export OpenBoundary, PeriodicBoundary, CollisionBoundary, MovingBoundary + +""" + North<:AbstractDirection + +A simple subtype of [`AbstractDirection`](@ref) used for parametrically typing a subtype of +[`AbstractBoundary`](@ref) if that boundary is the northern boundary in a rectangular domain. +""" +struct North<:AbstractDirection end + +#= + _boundary_info_from_extent(North, FT, x0, xf, y0, yf) + +Create a rectangular polygon representing the Northern boundary of a rectangular domain +defined by the x-extent (x0, xf) and the y-extent (y0, yf). The points of this boundary will +be of float type `FT`. If the length of the boundary x-extent is `Lx = xf - x0` then the +x-extent of the boundary polygon will range from `x0 - Lx/2` to `xf + Lx/2` in the x-direction. +If the length of the grid y-extent is `Ly = yf - y0` then the boundary polygon will range from +`yf` to `yf + Ly/2` in the y-direction. This will create overlap with other boundary walls, +if all created from the grid, making sure all floes connect with boundaries at edges of +the domain. + +Also return value `yf` as the value representing the edge of the boundary connecting with +the edge of the domain. +=# +function _boundary_info_from_extent(::Type{North}, ::Type{FT}, x0, xf, y0, yf) where FT + Δx, Δy = (xf - x0)/2, (yf - y0)/2 + poly = _make_bounding_box_polygon(FT, x0 - Δx, xf + Δx, yf, yf + Δy) + return poly, FT(yf) +end + +#= For floe interaction points (fpoints) that cross over the North boundary, zero-out all +forces in in direction not perpendicular to North boundary wall. =# +function _normal_direction_correct!(forces, fpoints, boundary::AbstractBoundary{North, FT}) where FT + forces[fpoints[:, 2] .>= boundary.val, 1] .= FT(0.0) + return +end + +""" + South<:AbstractDirection + +A simple subtype of [`AbstractDirection`](@ref) used for parametrically typing a subtype of +[`AbstractBoundary`](@ref) if that boundary is the southern boundary in a rectangular domain. +""" +struct South<:AbstractDirection end + +#= + _boundary_info_from_extent(South, FT, x0, xf, y0, yf) + +Create a rectangular polygon representing the Southern boundary of a rectangular domain +defined by the x-extent (x0, xf) and the y-extent (y0, yf). The points of this boundary will +be of float type `FT`. If the length of the grid x-extent is `Lx = xf - x0` then the +x-extent of the boundary polygon will range from `x0 - Lx/2` to `xf + Lx/2` in the x-direction. +If the length of the grid y-extent is `Ly = yf - y0` then the boundary polygon will range from +`y0 - Ly/2` to `y0` in the y-direction. This will create overlap with other boundary walls, +if all created from the grid, making sure all floes connect with boundaries at edges of +the domain. + +Also return value `y0` as the value representing the edge of the boundary connecting with +the edge of the domain. +=# +function _boundary_info_from_extent(::Type{South}, ::Type{FT}, x0, xf, y0, yf) where FT + Δx, Δy = (xf - x0)/2, (yf - y0)/2 + poly = _make_bounding_box_polygon(FT, x0 - Δx, xf + Δx, y0 - Δy, y0) + return poly, FT(y0) +end + +#= For floe interaction points (fpoints) that cross over the South boundary, zero-out all +forces in in direction not perpendicular to South boundary wall. =# +function _normal_direction_correct!(forces, fpoints, boundary::AbstractBoundary{South, FT}) where FT + forces[fpoints[:, 2] .<= boundary.val, 1] .= FT(0.0) + return +end + +""" + East<:AbstractDirection + + +A simple subtype of [`AbstractDirection`](@ref) used for parametrically typing a subtype of +[`AbstractBoundary`](@ref) if that boundary is the eastern boundary in a rectangular domain. +""" +struct East<:AbstractDirection end + +#= + _boundary_info_from_extent(East, FT, x0, xf, y0, yf) + +Create a rectangular polygon representing the Eastern boundary of a rectangular domain +defined by the x-extent (x0, xf) and the y-extent (y0, yf). The points of this boundary will +be of float type `FT`. If the length of the grid x-extent is `Lx = xf - x0` then the +x-extent of the boundary polygon will range from `xf` to `xf + Lx/2` in the x-direction. +If the length of the grid y-extent is `Ly = yf - y0` then the boundary polygon will range from +`y0 - Ly/2` to `yf + Ly/2` in the y-direction. This will create overlap with other boundary walls, +if all created from the grid, making sure all floes connect with boundaries at edges of +the domain. + +Also return value `xf` as the value representing the edge of the boundary connecting with +the edge of the domain. +=# +function _boundary_info_from_extent(::Type{East}, ::Type{FT}, x0, xf, y0, yf) where FT + Δx, Δy = (xf - x0)/2, (yf - y0)/2 + poly = _make_bounding_box_polygon(FT, xf, xf + Δx, y0 - Δy, yf + Δy) + return poly, FT(xf) +end + +#= For floe interaction points (fpoints) that cross over the East boundary, zero-out all +forces in in direction not perpendicular to East boundary wall. =# +function _normal_direction_correct!(forces, fpoints, boundary::AbstractBoundary{East, FT}) where FT + forces[fpoints[:, 1] .>= boundary.val, 2] .= FT(0.0) + return +end + +""" + West<:AbstractDirection + + +A simple subtype of [`AbstractDirection`](@ref) used for parametrically typing a subtype of +[`AbstractBoundary`](@ref) if that boundary is the western boundary in a rectangular domain. +""" +struct West<:AbstractDirection end + +#= + _boundary_info_from_extent(West, FT, x0, xf, y0, yf) + +Create a rectangular polygon representing the Western boundary of a rectangular domain +defined by the x-extent (x0, xf) and the y-extent (y0, yf). The points of this boundary will +be of float type `FT`. If the length of the grid x-extent is `Lx = xf - x0` then the +x-extent of the boundary polygon will range from `x0 - Lx/2` to `x0` in the x-direction. +If the length of the grid y-extent is `Ly = yf - y0` then the boundary polygon will range from +`y0 - Ly/2` to `yf + Ly/2` in the y-direction. This will create overlap with other boundary walls, +if all created from the grid, making sure all floes connect with boundaries at edges of +the domain. + +Also return value `x0` as the value representing the edge of the boundary connecting with +the edge of the domain. +=# +function _boundary_info_from_extent(::Type{West}, ::Type{FT}, x0, xf, y0, yf) where FT + Δx, Δy = (xf - x0)/2, (yf - y0)/2 + poly = _make_bounding_box_polygon(FT, x0 - Δx, x0, y0 - Δy, yf + Δy) + return poly, FT(x0) +end + +#= For floe interaction points (fpoints) that cross over the West boundary, zero-out all +forces in in direction not perpendicular to West boundary wall. =# +function _normal_direction_correct!(forces, fpoints, boundary::AbstractBoundary{West, FT}) where FT + forces[fpoints[:, 1] .<= boundary.val, 2] .= FT(0.0) + return +end + +# Concrete subtype of AbstractBoundary - see documentation below +struct OpenBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::StaticQuadrilateral{FT} + val::FT +end + +""" + OpenBoundary{D, FT} <: AbstractBoundary{D, FT} + +A concrete subtype of [`AbstractBoundary`](@ref) that removes a floe from the simulation if +any of the floe's vertices overlap with the `OpenBoundary`. This is meant to simulate the +floe floating out of the simulation. The floe is immeditaly removed as there might not be +ocean `u` and `v` values outside of the domian and thus we wouldn't know how to evolve the +floe's trajectory. This boundary is direction `D` ([`AbstractDirection`](@ref)) of the +domain and has field data of type `FT <: AbstractFloat`. + +## _Fields_ +- $BOUNDARY_POLY_DEF +- $VAL_DEF + +Here is how to construct an `OpenBoundary`: + + OpenBoundary(D, [FT = Float64]; grid = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing) + +The user must specify which [`AbstractDirection`](@ref) the boundary is. The user also has an opportunity to specify which +float type `Float64` or `Float32` will be used for field data. The user then must either provide +an grid object ([`AbstractRectilinearGrid`](@ref)) for the domain to align with the grid edges, or provide +four values (`x0`, `xf`, `y0`, and `yf`) that define the x and y-extents the user wants for the +domain. Note: if the user chooses to specify the domain extents, they still must be within the +grid in order to make a valid [`Model`](@ref). + +## _Positional arguments_ +- $D_DEF +- $FT_DEF + +## _Keyword arguments_ +- $GRID_DEF +- $X0_DEF +- $XF_DEF +- $Y0_DEF +- $YF_DEF + +_Note:_ the user must either provide a `grid` OR all four (4) of `x0`, `xf`, `y0`, and `yf`. + +## _Examples_ +- Defining a Northern `OpenBoundary` with Float64 (default type) data using the `grid` keyword. +```jldoctest +julia> g = RegRectilinearGrid(x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20); + +julia> OpenBoundary(North; grid = g) +OpenBoundary{North, Float64} + ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0) + ∟val is 500000.0 +``` +- Defining a Southern `OpenBoundary` with Float32 data using the `x0`, `xf`, `y0` and `yf` keywords. +```jldoctest +julia> OpenBoundary(South, Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5) +OpenBoundary{South, Float32} + ⊢polygon points are defined by the following set: (750000.0f0, 0.0f0), (750000.0f0, -250000.0f0), (-250000.0f0, 0.0f0), (-250000.0f0, -250000.0f0) + ∟val is 0.0 +``` +""" +function OpenBoundary(::Type{D}, ::Type{FT} = Float64; grid = nothing, + x0 = nothing, xf = nothing, y0 = nothing, yf = nothing, +) where {D, FT} + if !isnothing(grid) + x0, xf, y0, yf = _get_grid_extent(grid) + elseif isnothing(x0) || isnothing(xf) || isnothing(y0) || isnothing(yf) + throw(ArgumentError("To create an OpenBoundary, either provide a grid or a x0, xf, y0, AND yf.")) + end + poly, val = _boundary_info_from_extent(D, FT, x0, xf, y0, yf) + return OpenBoundary{D, FT}(poly, val) +end + +# Pretty printing for OpenBoundary showing key types and fields +function Base.show(io::IO, open_bound::OpenBoundary{D, FT}) where {D, FT} + overall_summary = "OpenBoundary{$D, $FT}" + points_summary, val_summary = show_boundary_poly_val_strings(open_bound) + print(io, overall_summary, "\n", + " ⊢", points_summary, "\n", + " ∟", val_summary) +end + +# Concrete subtype of AbstractBoundary - see documentation below +struct PeriodicBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::StaticQuadrilateral{FT} + val::FT +end + +""" + PeriodicBoundary <: AbstractBoundary + +A concrete subtype of [`AbstractBoundary`](@ref) that moves a floe from one side of the +domain to the opposite side of the domain ([`North`](@ref) to [`South`](@ref) and +[`East`](@ref) to [`West`](@ref) and visa versa) when its centroid crosses the +`PeriodicBoundary`, bringing the floe back into the domain. Due to this behavior, +`PeriodicBoundary` pairs are required to form a valid domain. This boundary is direction `D` +([`AbstractDirection`](@ref)) of the domain and has field data of type `FT <: AbstractFloat`. + +## _Fields_ +- $BOUNDARY_POLY_DEF +- $VAL_DEF + +Here is how to construct an `PeriodicBoundary`: + + PeriodicBoundary(D, [FT = Float64]; grid = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing) + +The user must specify which [`AbstractDirection`](@ref) the boundary is. The user also has an opportunity to specify which +float type `Float64` or `Float32` will be used for field data. The user then must either provide +an grid object ([`AbstractRectilinearGrid`](@ref)) for the domain to align with the grid edges, or provide +four values (`x0`, `xf`, `y0`, and `yf`) that define the x and y-extents the user wants for the +domain. Note: if the user chooses to specify the domain extents, they still must be within the +grid in order to make a valid [`Model`](@ref). + +## _Positional arguments_ +- $D_DEF +- $FT_DEF + +## _Keyword arguments_ +- $GRID_DEF +- $X0_DEF +- $XF_DEF +- $Y0_DEF +- $YF_DEF + +_Note:_ the user must either provide a `grid` OR all four (4) of `x0`, `xf`, `y0`, and `yf`. + +## _Examples_ +- Defining an Eastern `PeriodicBoundary` with Float32 data using the `grid` keyword. +```jldoctest +julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20); + +julia> PeriodicBoundary(East; grid = g) +PeriodicBoundary{East, Float64} + ⊢polygon points are defined by the following set: (750000.0, -250000.0), (500000.0, -250000.0), (750000.0, 750000.0), (500000.0, 750000.0) + ∟val is 500000.0 +``` +- Defining a Western `PeriodicBoundary` with Float64 data using the `x0`, `xf`, `y0` and `yf` keywords. +```jldoctest +julia> PeriodicBoundary(West, Float64; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5) +PeriodicBoundary{West, Float64} + ⊢polygon points are defined by the following set: (0.0, -250000.0), (-250000.0, 750000.0), (0.0, 750000.0), (-250000.0, -250000.0) + ∟val is 0.0 +``` +""" +function PeriodicBoundary(::Type{D}, ::Type{FT} = Float64; grid = nothing, + x0 = nothing, xf = nothing, y0 = nothing, yf = nothing, +) where {D, FT} + if !isnothing(grid) + x0, xf, y0, yf = _get_grid_extent(grid) + elseif isnothing(x0) || isnothing(xf) || isnothing(y0) || isnothing(yf) + throw(ArgumentError("To create an PeriodicBoundary, either provide a grid or an x0, xf, y0, AND yf.")) + end + poly, val = _boundary_info_from_extent(D, FT, x0, xf, y0, yf) + return PeriodicBoundary{D, FT}(poly, val) +end + +# Pretty printing for PeriodicBoundary showing key types and fields +function Base.show(io::IO, periodic_bound::PeriodicBoundary{D, FT}) where {D, FT} + overall_summary = "PeriodicBoundary{$D, $FT}" + points_summary, val_summary = show_boundary_poly_val_strings(periodic_bound) + print(io, overall_summary, "\n", + " ⊢", points_summary, "\n", + " ∟", val_summary) +end + +# Concrete type of AbstractBoundary - see documentation below +struct CollisionBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::StaticQuadrilateral{FT} + val::FT +end + +""" + CollisionBoundary <: AbstractBoundary + +A concrete subtype of [`AbstractBoundary`](@ref) that calculates collision forces of a floe +against the boundary if any of the floe's vertices overlap with the `CollisionBoundary`. +This is meant to simulate any barrier that might stop a floe from flowing into or out of a +given region, like the edges of a cove. With this type of wall, a `CollisionBoundary` is +treated as an immovable, unbreakable floe for the purposes of calculations. This boundary is +direction `D` ([`AbstractDirection`](@ref)) of the domain and has field data of type +`FT <: AbstractFloat`. + +## _Fields_ +- $BOUNDARY_POLY_DEF +- $VAL_DEF + +Here is how to construct an `CollisionBoundary`: + + CollisionBoundary(D, [FT = Float64]; grid = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing) + +The user must specify which [`AbstractDirection`](@ref) the boundary is. The user also has an opportunity to specify which +float type `Float64` or `Float32` will be used for field data. The user then must either provide +an grid object ([`AbstractRectilinearGrid`](@ref)) for the domain to align with the grid edges, or provide +four values (`x0`, `xf`, `y0`, and `yf`) that define the x and y-extents the user wants for the +domain. Note: if the user chooses to specify the domain extents, they still must be within the +grid in order to make a valid [`Model`](@ref). + +## _Positional arguments_ +- $D_DEF +- $FT_DEF + +## _Keyword arguments_ +- $GRID_DEF +- $X0_DEF +- $XF_DEF +- $Y0_DEF +- $YF_DEF + +_Note:_ the user must either provide a `grid` OR all four (4) of `x0`, `xf`, `y0`, and `yf`. + +## _Examples_ +- Defining an Northern `CollisionBoundary` with Float64 data using the `grid` keyword. +```jldoctest +julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20); + +julia> CollisionBoundary(North; grid = g) +CollisionBoundary{North, Float64} + ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0) + ∟val is 500000.0 +``` +- Defining a Western `CollisionBoundary` with Float64 data using the `x0`, `xf`, `y0` and `yf` keywords. +```jldoctest +julia> CollisionBoundary(West, Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5) +CollisionBoundary{West, Float32} + ⊢polygon points are defined by the following set: (0.0f0, -250000.0f0), (-250000.0f0, 750000.0f0), (0.0f0, 750000.0f0), (-250000.0f0, -250000.0f0) + ∟val is 0.0 +``` +""" +function CollisionBoundary(::Type{D}, ::Type{FT} = Float64; grid = nothing, + x0 = nothing, xf = nothing, y0 = nothing, yf = nothing, +) where {D, FT} + if !isnothing(grid) + x0, xf, y0, yf = _get_grid_extent(grid) + elseif isnothing(x0) || isnothing(xf) || isnothing(y0) || isnothing(yf) + throw(ArgumentError("To create an CollisionBoundary, either provide a grid or an x0, xf, y0, AND yf.")) + end + poly, val = _boundary_info_from_extent(D, FT, x0, xf, y0, yf) + return CollisionBoundary{D, FT}(poly, val) +end + +# Pretty printing for CollisionBoundary showing key types and fields +function Base.show(io::IO, collision_bound::CollisionBoundary{D, FT}) where {D, FT} + overall_summary = "CollisionBoundary{$D, $FT}" + points_summary, val_summary = show_boundary_poly_val_strings(collision_bound) + print(io, overall_summary, "\n", + " ⊢", points_summary, "\n", + " ∟", val_summary) +end + +# Concrete type of AbstractBoundary - see documentation below +mutable struct MovingBoundary{D, FT}<:AbstractBoundary{D, FT} + poly::StaticQuadrilateral{FT} + val::FT + u::FT + v::FT +end + +""" + MovingBoundary <: AbstractBoundary + +A concrete subtype of [`AbstractBoundary`](@ref) that can provide a compressive or shear +force on the floes within the domain by moving parallel or perpendicular to the domain. +The `MovingBoundary` calcuates the forces on the floe exactly the same as a +[`CollisionBoundary`](@ref), by acting as a immovable, unbreakable floe. This boundary is +direction `D` ([`AbstractDirection`](@ref)) of the domain and has field data of type +`FT <: AbstractFloat`. + +If a `North` or `South` wall has a non-zero `v` velocity this will provide a compressive +stress as the floe moves towards the center of the domain along the vector `x = cx` where +`(cx, cy)` is the centroid of the domain (or decompresses with opposite sign velocities). +The `poly` and `val` fields are updated to represent the movement of the domain at the +user-provided velocities, `u` and `v` m/s. + +Alternatively, if a `North` or `South` wall has a non-zero `u` velocity this will provide a +shear stress as the domain "moves" along the line `y = y0` (for `North`) or `y = yf` +(for South). In this case, this only changes the frictional forces and we do not need up +actually change the `poly` or `val` fields. + +## _Fields_ +- $BOUNDARY_POLY_DEF +- $VAL_DEF +- `u::FT`: boundary's u-velocity +- `v::FT`: boundary's v-velocity + +Here is how to construct an `MovingBoundary`: + + MovingBoundary(D, [FT = Float64]; u = 0.0, v = 0.0, grid = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing) + +The user must specify which [`AbstractDirection`](@ref) the boundary is. The user also has an opportunity to specify which +float type `Float64` or `Float32` will be used for field data. The user then must either provide +an grid object ([`AbstractRectilinearGrid`](@ref)) for the domain to align with the grid edges, or provide +four values (`x0`, `xf`, `y0`, and `yf`) that define the x and y-extents the user wants for the +domain. Note: if the user chooses to specify the domain extents, they still must be within the +grid in order to make a valid [`Model`](@ref). + +## _Positional arguments_ +- $D_DEF +- $FT_DEF + +## _Keyword arguments_ +- `u::FT`: boundary's u-velocity +- `v::FT`: boundary's v-velocity +- $GRID_DEF +- $X0_DEF +- $XF_DEF +- $Y0_DEF +- $YF_DEF + +_Note:_ the user must either provide a `grid` OR all four (4) of `x0`, `xf`, `y0`, and `yf`. +If the user does not provide values for `u` or `v`, the boundary will not move. + +## _Examples_ +- Defining an Northern `MovingBoundary` with Float64 data using the `grid` keyword. Assigning u-velocity of 0.5m/s. +```jldoctest +julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20); + +julia> MovingBoundary(North; u = 0.5, grid = g) +MovingBoundary{North, Float64} + ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0) + ⊢val is 500000.0 + ⊢u-velocity of 0.5 m/s + ∟v-velocity of 0.0 m/s +``` +- Defining a Southern `MovingBoundary` with Float32 data using the `x0`, `xf`, `y0` and `yf` keywords. +Assigning u-velocity of 0.3 m/s and v-velocity of 0.25 m/s +```jldoctest +julia> MovingBoundary(South, Float32; u = 0.3, v = 0.25, x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5) +MovingBoundary{South, Float32} + ⊢polygon points are defined by the following set: (750000.0f0, 0.0f0), (750000.0f0, -250000.0f0), (-250000.0f0, 0.0f0), (-250000.0f0, -250000.0f0) + ⊢val is 0.0 + ⊢u-velocity of 0.3 m/s + ∟v-velocity of 0.25 m/s +``` +""" +function MovingBoundary( + ::Type{D}, # Directions + ::Type{FT} = Float64; + u = 0.0, v = 0.0, grid = nothing, + x0 = nothing, xf = nothing, y0 = nothing, yf = nothing, +) where {D, FT} + if !isnothing(grid) + x0, xf, y0, yf = _get_grid_extent(grid) + elseif isnothing(x0) || isnothing(xf) || isnothing(y0) || isnothing(yf) + throw(ArgumentError("To create an MovingBoundary, either provide a grid or x0, xf, y0, AND yf.")) + end + if u == v == 0 + @warn "MovingBoundary velocities are both zero. Boundary will not move. Use keywords u and v to set velocities." + end + poly, val = _boundary_info_from_extent(D, FT, x0, xf, y0, yf) + return MovingBoundary{D, FT}(poly, val, u, v) +end + +# Return MovingBoundary velocity as assigned by the user. +_get_velocity(boundary::MovingBoundary, _, _) = (boundary.u, boundary.v) + +#= Move North/South MovingBoundary according to v-velocity and length of timestep by +changing the `poly` and `val` fields. =# +function _update_boundary!(boundary::MovingBoundary{D, FT}, Δt) where {D <: Union{North, South}, FT} + (x0, xf), (y0, yf) = GI.extent(boundary.poly) + Δy = boundary.v * Δt + y0, yf = y0 + Δy, yf + Δy + boundary.poly = _make_bounding_box_polygon(FT, x0, xf, y0, yf) + boundary.val += Δy + return +end + +#= Move East/West MovingBoundary according to u-velocity and length of timestep by +changing the `poly` and `val` fields. =# +function _update_boundary!(boundary::MovingBoundary{D, FT}, Δt) where {D <: Union{East, West}, FT} + (x0, xf), (y0, yf) = GI.extent(boundary.poly) + Δx = boundary.u * Δt + x0, xf = x0 + Δx, xf + Δx + boundary.poly = _make_bounding_box_polygon(FT, x0, xf, y0, yf) + boundary.val += Δx + return +end + +# Pretty printing for CollisionBoundary showing key types and fields +function Base.show(io::IO, moving_bound::MovingBoundary{D, FT}) where {D, FT} + overall_summary = "MovingBoundary{$D, $FT}" + points_summary, val_summary = show_boundary_poly_val_strings(moving_bound) + u_velocity_summary = "u-velocity of $(moving_bound.u) m/s" + v_velocity_summary = "v-velocity of $(moving_bound.v) m/s" + print(io, overall_summary, "\n", + " ⊢", points_summary, "\n", + " ⊢", val_summary, "\n", + " ⊢", u_velocity_summary, "\n", + " ∟", v_velocity_summary) +end + +# Union of all non-peridic boundary types to use as shorthand for dispatch. +const NonPeriodicBoundary = Union{OpenBoundary, CollisionBoundary, MovingBoundary} + +const StationaryBoundary{D, FT} = Union{OpenBoundary{D, FT}, CollisionBoundary{D, FT}, PeriodicBoundary{D, FT}} where {D, FT} + +# Return zero-velocities for StationaryBoundary. +_get_velocity(::StationaryBoundary{D, FT}, _, _) where {D, FT} = (zero(FT), zero(FT)) + +# Do NOT update a StationaryBoundary during timestep. +_update_boundary!(::StationaryBoundary, Δt) = return + +#= Boundaries across from one another (North and South OR East and West) are "periodically +compatible" if they are either both periodic, or if neither is periodic. This is because if +a floe exits a periodic boundary, it must be able to re-enter the opposite boundary to +fulfill its definition of periodic. This function is used to define valid domain walls.=# +_periodic_compat(::PeriodicBoundary, ::PeriodicBoundary) = true +_periodic_compat(::PeriodicBoundary, ::NonPeriodicBoundary) = false +_periodic_compat(::NonPeriodicBoundary, ::PeriodicBoundary) = false +_periodic_compat(::NonPeriodicBoundary, ::NonPeriodicBoundary) = true diff --git a/src/simulation_components/domain_components/domains.jl b/src/simulation_components/domain_components/domains.jl new file mode 100644 index 0000000..e9170ea --- /dev/null +++ b/src/simulation_components/domain_components/domains.jl @@ -0,0 +1,163 @@ +# Domain definition (combines 4 boundaries and topography) +export Domain + +struct Domain{FT, NB, SB, EB, WB, TT} + north::NB + south::SB + east::EB + west::WB + topography::TT + + function Domain{FT, NB, SB, EB, WB, TT}(north::NB, south::SB, east::EB, west::WB, topography::TT) where { + FT<:AbstractFloat, + NB<:AbstractBoundary{North, FT}, + SB<:AbstractBoundary{South, FT}, + EB<:AbstractBoundary{East, FT}, + WB<:AbstractBoundary{West, FT}, + TT<:TopographyField{FT}, + } + if !_periodic_compat(north, south) + throw(ArgumentError("North and south boundary walls are not \ + periodically compatable as only one of them is periodic.")) + elseif !_periodic_compat(east, west) + throw(ArgumentError("East and west boundary walls are not \ + periodically compatable as only one of them is periodic.")) + elseif north.val < south.val + throw(ArgumentError("North boundary value is less than south \ + boundary value.")) + elseif east.val < west.val + throw(ArgumentError("East boundary value is less than west \ + boundary value.")) + end + new{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) + end +end + +""" + Domain{FT, NB, SB, EB, WB, TT} + +A simulation `Domain` holds four (4) boundary elements (concrete subtypes of +[`AbstractBoundary`](@ref)) and a list (potentially empty) of [`TopographyElement`](@ref) +objects. There must be a wall typed by each cardinal direction, and each wall must hold data of the +same float type `FT`. Additionally, the [`North`](@ref) boundary should be "higher" (on the +y-axis) than the [`South`](@ref) boundary and the [`East`](@ref) should be futher to the +"right" (on the x-axis) than the [`West`](@ref) boundary. Aditionally, the boundary walls +should be overlapping as detailed in [`AbstractBoundary`](@ref). + +Additionally, the set of walls must be periodically compatible. This means that pairs of +opposite boundaries (`North` and `South` AND `East` and `West`) both need to be periodic if +one of them is periodic. This is because if a floe exits a periodic boundary, it must be +able to re-enter the opposite boundary to fulfill its definition of periodic. + +## _Fields_ +- `north::NB`: Northern boundary where `NB <: AbstractBoundary{North, FT}` +- `south::SB`: Southern boundary where `SB <: AbstractBoundary{South, FT}` +- `east::EB`: Eastern boundary where `EB <: AbstractBoundary{East, FT}` +- `west::WB`: Western boundary where `WB <: AbstractBoundary{West, FT}` +- `topography::tT`: Field of topography elements where `TT <: TopographyField{FT}` + +Notes: +- All `FT` values above must be the same float type to form a valid domain. +- The code depends on the boundaries forming a rectangle oriented along the +cartesian grid. Other shapes/orientations are not supported at this time. + +Here is how to construct an `MovingBoundary`: + + Domain(; north::NB, south::SB, east::EB, west::WB, topography::TT = nothing) + +The user must provide the four (4) boundaries. They then have an option to provide topography. +If no topography is provided, an empty `TopographyField` will be created with the same `FT` +as the first of the boundaries (which should be shared amoung boundaries). + +## _Keyword arguments_ +- `north::NB`: domain's northern boundary +- `south::SB`: domain's southern boundary +- `east::EB`: domain's eastern boundary +- `west::WB`: domain's western boundary +- `topography::TT`: domain's topography field, if there is one, else an empty field is created + +## _Examples_ +- Creating a `Domain` with _NO_ `topography` +```jldoctest domain +julia> grid = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20); + +julia> north = OpenBoundary(North, Float64; grid); + +julia> south = OpenBoundary(South, Float64; grid); + +julia> east = PeriodicBoundary(East, Float64; grid); + +julia> west = PeriodicBoundary(West, Float64; grid); + +julia> Domain(; north, south, east, west) +Domain + ⊢Northern boundary of type OpenBoundary{North, Float64} + ⊢Southern boundary of type OpenBoundary{South, Float64} + ⊢Eastern boundary of type PeriodicBoundary{East, Float64} + ⊢Western boundary of type PeriodicBoundary{West, Float64} + ∟0-element TopograpahyElement{Float64} list +``` + +- Creating a `Domain` with `topography` +```jldoctest domain +julia> import GeoInterface as GI; + +julia> topo_polys = [GI.Polygon([[(1e4, 1e4), (1e4, 3e4), (3e4, 1e4), (1e4, 1e4)]]), GI.Polygon([[(8e4, 8e4), (8e4, 9e4), (9e4, 9e4), (9e4, 8e4), (8e4, 8e4)]])]; + +julia> topography = initialize_topography_field(; polys = topo_polys); + +julia> Domain(; north, south, east, west, topography) +Domain + ⊢Northern boundary of type OpenBoundary{North, Float64} + ⊢Southern boundary of type OpenBoundary{South, Float64} + ⊢Eastern boundary of type PeriodicBoundary{East, Float64} + ⊢Western boundary of type PeriodicBoundary{West, Float64} + ∟2-element TopograpahyElement{Float64} list +``` +""" +function Domain(; + north::NB, south::SB, east::EB, west::WB, # need four boundary walls + topography::TT = nothing, # optional topography input +) where {NB, SB, EB, WB, TT} + FT = _get_boundary_float_type(north) + if isnothing(topography) # create an empty topography field + return Domain(; north, south, east, west, topography = initialize_topography_field(FT)) + end + return Domain{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) +end + +# Helper function to access `FT` type from boundary +_get_boundary_float_type(::AbstractBoundary{<:AbstractDirection, FT}) where FT = FT + +# Pretty printing for CollisionBoundary showing key types and fields +function Base.show(io::IO, domain::Domain{FT, NB, SB, EB, WB, TT}) where {FT, NB, SB, EB, WB, TT} + overall_summary = "Domain" + north_summary = "Northern boundary of type $NB" + south_summary = "Southern boundary of type $SB" + east_summary = "Eastern boundary of type $EB" + west_summary = "Western boundary of type $WB" + topo_summary = "$(length(domain.topography))-element TopograpahyElement{$FT} list" + + print(io, overall_summary, "\n", + " ⊢", north_summary, "\n", + " ⊢", south_summary, "\n", + " ⊢", east_summary, "\n", + " ⊢", west_summary, "\n", + " ∟", topo_summary) +end + +#= =# +function get_domain_element(domain, idx) + if idx == -1 + return domain.north + elseif idx == -2 + return domain.south + elseif idx == -3 + return domain.east + elseif idx == -4 + return domain.west + else + topo_idx = -(idx + 4) + return get_floe(domain.topography, topo_idx) + end +end diff --git a/src/simulation_components/domain_components/topography.jl b/src/simulation_components/domain_components/topography.jl new file mode 100644 index 0000000..09a5db5 --- /dev/null +++ b/src/simulation_components/domain_components/topography.jl @@ -0,0 +1,180 @@ +export TopographyElement # Topographic element within domain +export initialize_topography_field # Function to create topography field for user + +# Concrete subtype of AbstractDomainElement - see documentation below +struct TopographyElement{FT}<:AbstractDomainElement{FT} + poly::Polys{FT} + centroid::Vector{FT} + rmax::FT +end + +""" + TopographyElement{FT}<:AbstractDomainElement{FT} + +A concrete subtype of [`AbstractDomainElement`](@ref) that represent topography in a given +simulation. These `TopographyElement` objects will be treated as be treated as islands or +coastline within the model in that they will not move or break due to floe interactions, but +they will affect floes that collide with them. + +Similar to [`Floe`](@ref) objects, we store the shape of each `TopographyElement` with a +`poly` field. We also store the `centroid` and the shape's maximum radius (`rmax`) in order +to calculate simple interaction checks with Floes. + +## _Fields_ +- $POLY_DEF +- $CENTROID_DEF +- $RMAX_DEF + +Here is how to construct an `TopographyElement`: + + TopographyElement([FT = Float64]; poly::Polys) + +The user must provide a `GeoInterface` polygon of the specific type defined by `Polys`. +The user also has an opportunity to specify which float type `Float64` or `Float32` will be +used for field data. + +## _Positional arguments_ +- $FT_DEF + +## _Keyword arguments_ +- $POLY_DEF + +## _Note_: The user should **NOT** be using this constructor. The user should create a topography +field using [`initialize_topography_field`](@ref) rather than creating individual topography +elements. This allows the abstraction of the use of the `StructArrays` package away from the +user. + +## _Examples_ +```jldoctest topography +julia> import GeoInterface as GI; + +julia> poly = GI.Polygon([[(0.0, 0.0), (0.0, 1e3), (1e3, 1e3), (1e3, 0.0), (0.0, 0.0)]]); + +julia> TopographyElement(Float64; poly) +TopographyElement{Float64} + ⊢centroid is (500.0, 500.0) in meters + ∟maximum radius is 707.1067811865476 meters +``` + +```jldoctest topography +julia> TopographyElement(Float32; poly) +TopographyElement{Float32} + ⊢centroid is (500.0, 500.0) in meters + ∟maximum radius is 707.1068 meters +``` +""" +function TopographyElement(::Type{FT} = Float64; poly::Polys) where FT + # Clean up polygon and calculate centroid and maximum radius + poly = GO.ClosedRing()(poly) + rmholes!(poly) + centroid = collect(GO.centroid(poly, FT)) # TODO: Remove collect once type is changed + rmax = calc_max_radius(poly, centroid, FT) + return TopographyElement{FT}(poly, centroid, rmax) +end + +# Return zero-velocities for topography elements +_get_velocity(::TopographyElement{FT}, _, _) where {FT} = (zero(FT), zero(FT)) + +# No forces should be zero-ed out in collidions with topography elements, unlike boundaries. +_normal_direction_correct!(_, _, ::TopographyElement) = return + +# Pretty printing for TopographyElement showing key types and fields +function Base.show(io::IO, topo_element::TopographyElement{FT}) where FT + overall_summary = "TopographyElement{$FT}" + centroid_summary = "centroid is ($(GI.x(topo_element.centroid)), $(GI.y(topo_element.centroid))) in meters" + rmax_summary = "maximum radius is $(topo_element.rmax) meters" + print(io, overall_summary, "\n", + " ⊢", centroid_summary, "\n", + " ∟", rmax_summary) +end + +""" + initialize_topography_field([::Type{FT} = Float64]; polys = nothing, coords = nothing) + +This function allows for easy initialization of a field of [`TopographyElement(@ref)s and +collects them into a `StructArray` so that they can be passed into a [`Model`](@ref). + +This is the suggested way to create a topography field for a simulation. Do NOT construct +individual `TopographyElement` objects as that method does not correct the field as a whole +to ensure no topography polygons overlap and does not create a struct array that can be +passed to a `Model`. + +The user can create a topography field by either passing a list of polygons or by passing +a list of coordiantes, which will then be made into polygons. + +## _Positional arguments_ +- $FT_DEF + +## _Keyword arguments_ +- $POLY_LIST_DEF +- $COORDS_LIST_DEF + +_Note_: Topography field elements must not be intersecting, so the corrections within this +function may lead to more polygons than input to make a valid topography field. + +## _Examples_ +- Defining a topography field with coordinates +```jldoctest topo_field +julia> coords = [[[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)]], [[(10.0, 0.0), (10.0, 1.0), (11.0, 1.0), (10.0, 0.0)]]]; + +julia> initialize_topography_field(Float64; coords) +2-element TopographyField{Float64} list: + TopographyElement{Float64} + ⊢centroid is (0.3333333333333333, 0.6666666666666666) in meters + ∟maximum radius is 0.74535599249993 meters + TopographyElement{Float64} + ⊢centroid is (10.333333333333334, 0.6666666666666666) in meters + ∟maximum radius is 0.7453559924999301 meters +``` + +- Defining a topography field with polygons +```jldoctest topo_field +julia> import GeoInterface as GI; + +julia> polys = [GI.Polygon(c) for c in coords]; + +julia> initialize_topography_field(Float32; polys) +2-element TopographyField{Float32} list: + TopographyElement{Float32} + ⊢centroid is (0.33333334, 0.6666667) in meters + ∟maximum radius is 0.745356 meters + TopographyElement{Float32} + ⊢centroid is (10.333333, 0.6666667) in meters + ∟maximum radius is 0.74535626 meters +``` + +- Creating an empty topography field without polys or coords +```jldoctest +julia> initialize_topography_field(Float64) +0-element TopographyField{Float64} list +``` +""" +function initialize_topography_field(::Type{FT} = Float64; polys = nothing, coords = nothing) where FT + # Make sure input given (if given) is turned into a list of polygons + if isnothing(polys) && !isnothing(coords) + polys = [make_polygon(c) for c in coords] + elseif isnothing(polys) # & isnothing(coords) + return StructArray{TopographyElement{FT}}(undef, 0) + end + # Make sure list of polygons is non-intersecting --> these polygons will form topography + topo_multipoly = GO.DiffIntersectingPolygons()(GI.MultiPolygon(polys)) + # Make each polygon from multipolygon into a topography element + topo_field = StructArray{TopographyElement{FT}}(undef, GI.npolygon(topo_multipoly)) + for (i, poly) in enumerate(GI.getpolygon(topo_multipoly)) + topo_field[i] = TopographyElement(FT; poly) + end + return topo_field +end + +""" + TopographyField{FT} + +Alias for StructArray type with TopographyElement elements with data of type `FT`. +Note: this is not designed to be used by users, but is useful for dispatch by developers. +""" +const TopographyField{FT} = StructArray{TopographyElement{FT}} where FT + +# Pretty printing for TopographyField showing key types and elements +function Base.showarg(io::IO, ::TopographyField{FT}, toplevel) where FT + print(io, "TopographyField{$FT} list") +end diff --git a/src/simulation_components/grids.jl b/src/simulation_components/grids.jl new file mode 100644 index 0000000..e9ef443 --- /dev/null +++ b/src/simulation_components/grids.jl @@ -0,0 +1,219 @@ +export AbstractRectilinearGrid, RegRectilinearGrid + +# See CellFloes documentation below +struct CellFloes{FT<:AbstractFloat} + floeidx::Vector{Int} + Δx::Vector{FT} + Δy::Vector{FT} +end + +""" + CellFloes([FT = Float64]; floeidx, Δx, Δy) + +Constructor for struct that represents a single grid cell and accumulates the indices of +floes with area in that grid cell. This is used for two-way coupling to accumulate the +forces of floes in a grid cell on the ocean below it. Due to the prevalence of periodic +boundaries, the grid cell represented by a `CellFloes` object are centered on grid points +(translated by `Δx/2` and `Δy/2` in the x and y directions for a [`RegRectilinearGrid`](@ref)), +rather than on the grid cells defined by the grid object itself. Floes are recorded with +their index in the list of floes. Furthermore, in a model with periodic boundaries, some +floes may be in multiple grid cells on different edges of the domain if they pass through a +periodic boundary. In these cases, the floe "linked" to with its index must be translated by +a vector to get its "ghost" on the other side of the domain. This mon-integer translation +data is of float type `FT`. + +**Note**: If no keyword arguments are provide by the user, an `CellFloes` object with empty +fields will be created. This is the **standard useage** of these objects and they are added to +during the coupling step. If keyword arguments are provided, then all three must be provided +and each vector must be the same size. + +## _Positional arguments_ +- $FT_DEF +## _Keyword arguments_ +- `floeidx::Vector{Int}`: vector of floe indicies in the list of model floes for floes with area in the grid +- `Δx::Vector{FT}`: vector of x-translations for a floe at the corresponding index of the +`floeidx` vector to be in the cell represented with the `CellFloes` object. +- `Δy::Vector{FT}`: vector of y-translations for a floe at the corresponding index of the +`floeidx` vector to be in the cell represented with the `CellFloes` object. + +## _CellFloes Fields_ +- `floeidx`: see keyword arguments +- `Δx`: see keyword arguments +- `Δy`: see keyword arguments +""" +function CellFloes(::Type{FT} = Float64; floeidx = nothing, Δx = nothing, Δy = nothing) where {FT} + if isnothing(floeidx) || isnothing(Δx) || isnothing(Δy) + floeidx = Vector{Int}() + Δx = Vector{FT}() + Δy = Vector{FT}() + else + @assert length(floeidx) == length(Δx) == length(Δy) "A CellFloes object requires that all fields be vectors of the same length." + end + return CellFloes{FT}(floeidx, Δx, Δy) +end + +""" + empty!(cell::CellFloes) + +Empties each of the three vectors (`floeidx`, `Δx`, and `Δy`) within a [`CellFloes`](@ref) object. +""" +function Base.empty!(cell::CellFloes) + empty!(cell.floeidx) + empty!(cell.Δx) + empty!(cell.Δy) + return +end + +const GRID_DEF = "`grid::AbstractRectilinearGrid`: subtype of AbstractRectilinearGrid representing the grid used within a simulation" + +""" + YourGridType{FT} <: AbstractRectilinearGrid{FT} + +Each simulation run with Subzero.jl must be run on a grid. A grid defines the points at +which the ocean and atmosphere hold velocity data, as well as various other tracers, at grid +points. The user must choose an existing subtype of `AbstractRectilinearGrid` or implement a new +subtype to create a simulation. + +Each grid implementation must define the number and dimensions of each grid cell. Right now, +this assumes that the ocean and the atmosphere are on the same grid and that the grid is +rectangular. We might not want this to be true in the future. Furthermore, each concrete +implementation of `AbstractRectilinearGrid` that will be used to two-way couple with an ocean +or an atmosphere must have a field called `floe_locations` that is a matrix of +[`CellFloes`](@ref), with one element for each grid point. The user should not worry about +`CellFloes`, this is up to the developer to make sure that their grid contains and populates +this field. Again, in the future, this might be related to the ocean and/or atmosphere +rather than the grid. For more information, or if you're interested in working on this, see +issue [#107](@ref). + +## _API_ +The following methods must be implemented for all subtypes: +- `_get_grid_extent(grid::AbstractRectilinearGrid)`: + +The `_get_grid_extent` function takes in a concrete subtype of `AbstractRectilinearGrid` and +returns the grid's minimum x-value, maximum x-value, minimum y-value, and maximum y-value. + +Given the original code was written for [`RegRectilinearGrid`](@ref) objects, the dispatch +isn't fully implemented and other functions should be added to this list. If you were +interested in adding a new subtype of `AbstractRectilinearGrid`, see issue [#108](@ref). +""" +abstract type AbstractRectilinearGrid{FT<:AbstractFloat} end + +# Concrete subtype of AbstractRectilinearGrid - see below for documentation +struct RegRectilinearGrid{FT}<:AbstractRectilinearGrid{FT} + Nx::Int + Ny::Int + Δx::FT + Δy::FT + x0::FT + xf::FT + y0::FT + yf::FT + floe_locations::Matrix{CellFloes{FT}} +end + +""" + RegRectilinearGrid{FT} <: AbstractRectilinearGrid{FT} + +A concrete implementation of an [`AbstractRectilinearGrid`](@ref) that represents a tessellation of 2-dimensional +Euclidean space into `n`-by-`m` congruent rectangles. Fields that hold float data are of +type `FT`, a concrete subtype of `AbstractFloat`. + +- `Nx::Int`: number of grid cells in the x-direction +- `Ny::Int`: number of grid cells in the y-direction +- `Δx::FT`: grid cell width +- `Δy::FT`: grid cell height +- `x0::FT`: value of first x grid line +- `xf::FT`: value of final x grid line +- `y0::FT`: value of first y grid line +- `yf::FT`: value of final y grid line +- `floe_locations::Matrix{CellFloes}`: `Nx + 1` by `Ny + 1` matrix of [`CellFloes`](@ref) + +Here is how to construct a `RegRectilinearGrid`: + + RegRectilinearGrid([FT = Float64]; x0, xf, y0, yf, Nx = nothing, Ny = nothing, Δx = nothing, Δy = nothing) + +## _Positional arguments_ +- $FT_DEF +## _Keyword arguments_ +- `x0::FT`: value of first x grid line +- `xf::FT`: value of final x grid line +- `y0::FT`: value of first y grid line +- `yf::FT`: value of final y grid line +- `Nx::Int`: number of grid cells in the x-direction +- `Ny::Int`: number of grid cells in the y-direction +- `Δx::FT`: grid cell width +- `Δy::FT`: grid cell height + +_Note:_ the user MUST provide `x0`, `xf`, `y0`, and `yf`; the user then has the choice to +provide `Nx` and `Ny` OR `Δx` and `Δy`. If provided `Δx` doesn't evenly divide length +`xf-x0` or `Δy` doesn't evenly divide `yf-y0`, you won't get full size grid. The grid will +be "trimmed" to the nearest full grid square in both directions. + +## _Examples_ + +- Defining a `RegRectilinearGrid` using `Nx` and `Ny`. +```jldoctest +julia> grid = RegRectilinearGrid(; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 10) +RegRectilinearGrid{Float64} + ⊢x extent (0.0 to 500000.0) with 20 grid cells of size 25000.0 m + ∟y extent (0.0 to 500000.0) with 10 grid cells of size 50000.0 m +``` +- Defining a `RegRectilinearGrid` using `Δx` and `Δy`. +```jldoctest +julia> grid = RegRectilinearGrid(Float32; x0 = -5e5, xf = 5e5, y0 = 0.0, yf = 5e5, Δx = 5e4, Δy = 5e4) +RegRectilinearGrid{Float32} + ⊢x extent (-500000.0 to 500000.0) with 20 grid cells of size 50000.0 m + ∟y extent (0.0 to 500000.0) with 10 grid cells of size 50000.0 m +``` +- Error due to attemping to define a `RegRectilinearGrid` using `Δx` and `Ny`. +```jldoctest +julia> grid = RegRectilinearGrid(; x0 = -5e5, xf = 5e5, y0 = 0.0, yf = 5e5, Δx = 5e4, Ny = 10) +ERROR: ArgumentError: To create a RegRectilinearGrid, either provide Δx and Δy OR Nx and Ny. +[...] +``` +""" +function RegRectilinearGrid(::Type{FT} = Float64; + x0, xf, y0, yf, + Nx = nothing, Ny = nothing, + Δx = nothing, Δy = nothing, +) where FT <: AbstractFloat + # ensure x and y endpoints are properly ordered + if x0 > xf + @warn "$x0 can't be larger than $xf. Swapping the order of the grid x-endpoints" + xf, x0 = x0, xf + end + if y0 > yf + @warn "$y0 can't be larger than $yf. Swapping the order of the grid y-endpoints" + yf, y0 = y0, yf + end + # determine number and size of gridpoints based on provided arguments + if !isnothing(Nx) && !isnothing(Ny) # determine size of grid cells based on number + Δx = (xf - x0) / Nx + Δy = (yf - y0) / Ny + elseif !isnothing(Δx) && !isnothing(Δy) # determine number of grid cells based on size + Nx = floor(Int, (xf - x0) / Δx) + Ny = floor(Int, (yf - y0) / Δy) + # note: round grid size if requested cell size doesn't go evenly into endpoints + xf = x0 + Nx * Δx + yf = y0 + Ny * Δy + else + throw(ArgumentError("To create a RegRectilinearGrid, either provide Δx and Δy OR Nx and Ny.")) + end + # create cell floe list + floe_locations = [CellFloes(FT) for _ in 1:(Nx + 1), _ in 1:(Ny + 1)] + # create RegRectilinearGrid object + return RegRectilinearGrid{FT}(Nx, Ny, Δx, Δy, x0, xf, y0, yf, floe_locations) +end + +# Return the grid's minimum x-value, maximum x-value, minimum y-value, and maximum y-value +_get_grid_extent(grid::RegRectilinearGrid) = (grid.x0, grid.xf, grid.y0, grid.yf) + +# Pretty printing for RegRectilinearGrid showing key dimensions +function Base.show(io::IO, grid::RegRectilinearGrid{FT}) where FT + overall_summary = "RegRectilinearGrid{$FT}" + x_summary = "x extent ($(grid.x0) to $(grid.xf)) with $(grid.Nx) grid cells of size $(grid.Δx) m" + y_summary = "y extent ($(grid.y0) to $(grid.yf)) with $(grid.Ny) grid cells of size $(grid.Δy) m" + print(io, overall_summary, "\n", + " ⊢", x_summary, "\n", + " ∟", y_summary) +end diff --git a/src/simulation_components/model.jl b/src/simulation_components/model.jl index 6415775..cb267ea 100644 --- a/src/simulation_components/model.jl +++ b/src/simulation_components/model.jl @@ -9,11 +9,11 @@ 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 + grid Outputs: true if domain is within grid bounds, else false """ -function domain_in_grid(domain::Domain, grid::AbstractGrid) +function domain_in_grid(domain::Domain, grid::AbstractRectilinearGrid) northval = domain.north.val southval = domain.south.val eastval = domain.east.val @@ -46,7 +46,7 @@ run. """ struct Model{ FT<:AbstractFloat, # float type - GT<:AbstractGrid{FT}, # grid type + GT<:AbstractRectilinearGrid{FT}, # grid type DT<:Domain{ # domain type FT, <:AbstractBoundary, @@ -70,7 +70,7 @@ struct Model{ floes::FLT, ) where { FT<:AbstractFloat, - GT<:AbstractGrid{FT}, + GT<:AbstractRectilinearGrid{FT}, DT<:Domain{ FT, <:AbstractBoundary, @@ -106,7 +106,7 @@ struct Model{ floes::FLT, ) where { FT<:AbstractFloat, - GT<:AbstractGrid{FT}, + GT<:AbstractRectilinearGrid{FT}, DT<:Domain{ FT, <:AbstractBoundary, diff --git a/src/simulation_components/ocean_and_atmos.jl b/src/simulation_components/ocean_and_atmos.jl index 267aa47..195f862 100644 --- a/src/simulation_components/ocean_and_atmos.jl +++ b/src/simulation_components/ocean_and_atmos.jl @@ -151,7 +151,7 @@ end Construct model's ocean. Inputs: - grid model grid + 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 @@ -159,7 +159,7 @@ Output: Ocean with constant velocity and temperature on each grid line. """ Ocean{FT}( - grid::AbstractGrid, + grid::AbstractRectilinearGrid, u, v, temp, @@ -214,7 +214,7 @@ Atmos(args...) = Atmos{Float64}(args...) Construct model atmosphere of type FT. Inputs: - grid model's grid + 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 @@ -222,7 +222,7 @@ Output: Atmosphere of type FT with constant velocity and temperature over domain. """ Atmos{FT}( - grid::AbstractGrid, + grid::AbstractRectilinearGrid, u, v, temp, diff --git a/src/simulation_components/simulation.jl b/src/simulation_components/simulation.jl index 2396432..300a41e 100644 --- a/src/simulation_components/simulation.jl +++ b/src/simulation_components/simulation.jl @@ -48,7 +48,7 @@ The user can also define settings for each physical process. """ @kwdef struct Simulation{ FT<:AbstractFloat, - MT<:Model{FT, <:AbstractGrid, <:Domain}, + MT<:Model{FT, <:AbstractRectilinearGrid, <:Domain}, CT<:AbstractFractureCriteria, PT<:AbstractSubFloePointsGenerator{FT}, ST<:AbstractStressCalculator{FT}, diff --git a/test/qualitative_behavior.jl b/test/qualitative_behavior.jl index 51a63fa..8183d7f 100644 --- a/test/qualitative_behavior.jl +++ b/test/qualitative_behavior.jl @@ -18,12 +18,7 @@ const coarse_nx = 10 const coarse_ny = 10 # Setup for Simulations -grid = RegRectilinearGrid( - (-2.5e4, Lx), - (-2.5e4, Ly), - Δgrid, - Δgrid, -) +grid = RegRectilinearGrid(; x0 = -2.5e4, xf = Lx, y0 = -2.5e4, yf = Ly, Δx = Δgrid, Δy = Δgrid) zero_ocn = Ocean(grid, 0.0, 0.0, 0.0) meridional_ocn = Ocean(grid, 0.0, 1.0, 0.0) @@ -31,27 +26,19 @@ meridional_ocn = Ocean(grid, 0.0, 1.0, 0.0) zero_atmos = Atmos(grid, 0.0, 0.0, 0.0) zonal_atmos = Atmos(grid, -15.0, 0.0, 0.0) -open_domain_no_topo = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), +open_domain_no_topo = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) -topography = TopographyElement( - [[ - [2e4, 0.0], - [2e4, 2e4], - [2.5e4, 2e4], - [2.5e4, 0.0], - [2e4, 0.0], - ]]) -collision_domain_topo = Subzero.Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), - StructArray([topography]), +collision_domain_topo = Subzero.Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = CollisionBoundary(East; grid), + west = CollisionBoundary(West; grid), + topography = initialize_topography_field(; coords = [[[[2e4, 0.0], [2e4, 2e4], [2.5e4, 2e4], [2.5e4, 0.0], [2e4, 0.0]]]]) ) stationary_rect_floe = StructArray([Floe( @@ -218,20 +205,12 @@ Expected Behavior: populating a ghost floe which hits the topography element before bounding back through the western wall. """ -periodic_bounds_topo = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), - StructArray([TopographyElement( - [[ - [-1.5e4, 4.5e4], - [-1.5e4, 6.5e4], - [2.5e4, 6.5e4], - [2.5e4, 4.5e4], - [-1.5e4, 4.5e4], - ]], - )]), +periodic_bounds_topo = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), + topography = StructArray([TopographyElement(; coords = [[[-1.5e4, 4.5e4], [-1.5e4, 6.5e4], [2.5e4, 6.5e4], [2.5e4, 4.5e4], [-1.5e4, 4.5e4]]])]), ) p1_coords = [[ diff --git a/test/runtests.jl b/test/runtests.jl index 3212803..f4a59a2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,10 @@ import GeometryOps.GeoInterface as GI using Test @testset "Subzero.jl" begin + include("test_simulation_components/grids.jl") + include("test_simulation_components/domain_components/boundaries.jl") + include("test_simulation_components/domain_components/domains.jl") + include("test_simulation_components/domain_components/topography.jl") include("test_simulation_components/test_stress_calculators.jl") include("test_physical_processes/test_update_floe.jl") include("test_physical_processes/test_collisions.jl") diff --git a/test/test_conservation.jl b/test/test_conservation.jl index 001256f..f9c90b1 100644 --- a/test/test_conservation.jl +++ b/test/test_conservation.jl @@ -59,33 +59,26 @@ end @testset "Conservation of Energy and Momentum" begin Δt = 10 FT = Float64 - grid = RegRectilinearGrid( - (-2e4, 1e5), - (0, 1e5), - 1e4, - 1e4, - ) - collision_domain = Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), - ) - open_domain = Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), - ) - topo = TopographyElement( - [[[-1e4, 0.0], [-2e4, 1e4], [-1e4, 1e4], [-1e4, 0.0]]], - ) - open_domain_w_topography = Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), - StructVector([topo]) + grid = RegRectilinearGrid(; x0 = -2e4, xf = 1e5, y0 = 0.0, yf = 1e5, Δx = 1e4, Δy = 1e4) + collision_domain = Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = CollisionBoundary(East; grid), + west = CollisionBoundary(West; grid), + ) + open_domain = Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), + ) + topography = initialize_topography_field(; coords = [[[[-1e4, 0.0], [-2e4, 1e4], [-1e4, 1e4], [-1e4, 0.0]]]]) + open_domain_w_topography = Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), + topography ) rng = Xoshiro(1) floe1 = [[[2e4, 2e4], [2e4, 5e4], [5e4, 5e4], [5e4, 2e4], [2e4, 2e4]]] diff --git a/test/test_floe.jl b/test/test_floe.jl index 795ca42..36b7e9a 100644 --- a/test/test_floe.jl +++ b/test/test_floe.jl @@ -71,20 +71,20 @@ @test !any(Subzero.hashole.(floe_arr.coords)) # Test initialize_floe_field from coord list - grid = RegRectilinearGrid((-Lx, Lx), (-Ly, Ly), Δgrid, Δgrid) - nbound = CollisionBoundary(North, grid) - sbound = CollisionBoundary(South, grid) - ebound = CollisionBoundary(East, grid) - wbound = CollisionBoundary(West, grid) - domain_no_topo = Domain(nbound, sbound, ebound, wbound) + grid = RegRectilinearGrid(; x0 = -Lx, xf = Lx, y0 = -Ly, yf = Ly, Δx = Δgrid, Δy = Δgrid) + nbound = CollisionBoundary(North; grid) + sbound = CollisionBoundary(South; grid) + ebound = CollisionBoundary(East; grid) + wbound = CollisionBoundary(West; grid) + domain_no_topo = Domain(; north = nbound, south = sbound, east = ebound, west = wbound) island = [[[6e4, 4e4], [6e4, 4.5e4], [6.5e4, 4.5e4], [6.5e4, 4e4], [6e4, 4e4]]] topo1 = [[[-8e4, -8e4], [-8e4, 8e4], [-6e4, 8e4], [-5e4, 4e4], [-6e4, -8e4], [-8e4, -8e4]]] - domain_with_topo = Domain( - nbound, - sbound, - ebound, - wbound, - StructVector([TopographyElement(t) for t in [island, topo1]]) + domain_with_topo = Domain(; + north = nbound, + south = sbound, + east = ebound, + west = wbound, + topography = initialize_topography_field(; coords = [island, topo1]) ) topo_polys = Subzero.make_multipolygon([island, topo1]) @@ -102,12 +102,12 @@ @test all(floe_arr.id .== range(1, nfloes)) # From file with small domain -> floes outside of domain - small_grid = RegRectilinearGrid((-Lx/2, Lx/2), (-Ly/2, Ly), Δgrid, Δgrid) - small_domain_no_topo = Domain( - CollisionBoundary(North, small_grid), - CollisionBoundary(South, small_grid), - CollisionBoundary(East, small_grid), - CollisionBoundary(West, small_grid), + small_grid = RegRectilinearGrid(; x0 = -Lx/2, xf = Lx/2, y0 = -Ly/2, yf = Ly, Δx = Δgrid, Δy = Δgrid) + small_domain_no_topo = Domain(; + north = CollisionBoundary(North; grid = small_grid), + south = CollisionBoundary(South; grid = small_grid), + east = CollisionBoundary(East; grid = small_grid), + west = CollisionBoundary(West; grid = small_grid), ) floe_arr = (@test_logs (:warn, "Some floe centroids are out of the domain.") initialize_floe_field( FT, diff --git a/test/test_model.jl b/test/test_model.jl index befa77a..0aaa53b 100644 --- a/test/test_model.jl +++ b/test/test_model.jl @@ -1,97 +1,9 @@ -@testset "Model Creation" begin - # Grid Creation - @testset "Grid" begin - # Default constructor fails for non-matching dimensions - @test_throws ArgumentError Subzero.RegRectilinearGrid( - 80, - 50, - 0, - 1e5, - 0, - 1e5, - 1e3, # 1e5/80 = 1.25e3 ≂̸ 1e3 - 2e3, - [CellFloes{Float64}() for i in 1:81, j in 1:51], - ) - @test_throws ArgumentError Subzero.RegRectilinearGrid( - 80, - 50, - 0, - 1e5, - 0, - 1e5, - 1.25e3, - 1e3, # 1e5/50 = 2e3 ≂̸ 1e3 - [CellFloes{Float64}() for i in 1:81, j in 1:51], - ) - @test_throws ArgumentError Subzero.RegRectilinearGrid( - 80, - 50, - 0, - 1e5, - 0, - 1e5, - 1.25e3, - 2e3, - [CellFloes{Float64}() for i in 1:71, j in 1:51], # wrong dims - ) - - # Non-square grid using constructor with Δx and Δy - g1 = Subzero.RegRectilinearGrid( - (-10, 10), - (-8, 8), - 2, - 4, - ) - @test g1.Nx == 10 - @test g1.Ny == 4 - @test g1.x0 == -10 - @test g1.xf == 10 - @test g1.y0 == -8 - @test g1.yf == 8 - @test size(g1.floe_locations) == (11, 5) - @test typeof(g1) == Subzero.RegRectilinearGrid{Float64} - - # Non-square grid using constructor with Nx and Ny - g2 = Subzero.RegRectilinearGrid( - 10, - 4, - (-10, 10), - (-8, 8), - ) - @test g2.x0 == -10 - @test g2.xf == 10 - @test g2.y0 == -8 - @test g2.yf == 8 - @test g2.Δx == 2 - @test g2.Δy == 4 - @test size(g2.floe_locations) == (11, 5) - @test typeof(g2) == Subzero.RegRectilinearGrid{Float64} - - # Custom constructor Float32 and Float64 - @test typeof(Subzero.RegRectilinearGrid( - Float32, - (0, 10), - (0, 8), - 2, - 2, - )) == Subzero.RegRectilinearGrid{Float32} - @test typeof(Subzero.RegRectilinearGrid( - Float64, - (0, 10), - (0, 8), - 2, - 2, - )) == Subzero.RegRectilinearGrid{Float64} - end +using Test, Subzero +import StaticArrays as SA +@testset "Model Creation" begin @testset "Ocean" begin - g = Subzero.RegRectilinearGrid( - (0, 4e5), - (0, 3e5), - 1e4, - 1e4, - ) + g = Subzero.RegRectilinearGrid(; x0 = 0, xf = 4e5, y0 = 0, yf = 3e5, Δx = 1e4, Δy = 1e4) # Large ocean default constructor uocn = fill(3.0, g.Nx + 1, g.Ny + 1) vocn = fill(4.0, g.Nx + 1, g.Ny + 1) @@ -137,12 +49,8 @@ end @testset "Atmos" begin - g = Subzero.RegRectilinearGrid( - (0, 4e5), - (0, 3e5), - 1e4, - 1e4, - ) + g = Subzero.RegRectilinearGrid(; x0 = 0, xf = 4e5, y0 = 0, yf = 3e5, Δx = 1e4, Δy = 1e4) + # Large Atmos default constructor uatmos = fill(3.0, g.Nx + 1, g.Ny + 1) vatmos = fill(4.0, g.Nx + 1, g.Ny + 1) @@ -163,177 +71,6 @@ Subzero.Atmos{Float64} end - @testset "Boundaries" begin - # Boundaries using BoundaryCoords - FT = Float64 - g = Subzero.RegRectilinearGrid( - (0, 4e5), - (0, 3e5), - 1e4, - 1e4, - ) - b1 = Subzero.PeriodicBoundary(North, g) - b2 = Subzero.OpenBoundary(East, g) - b3 = Subzero.CollisionBoundary(West, g) - b4 = Subzero.PeriodicBoundary(South, g) - b5 = Subzero.MovingBoundary(South, g, 1.0, 2.0) - @test b1.val == 3e5 - @test typeof(b1) == Subzero.PeriodicBoundary{North, Float64} - @test b1.coords == [[[-2e5, 3e5], [-2e5, 4.5e5], [6e5, 4.5e5], [6e5, 3e5], [-2e5, 3e5]]] - @test b2.val == 4e5 - @test typeof(b2) == Subzero.OpenBoundary{East, Float64} - @test b2.coords == [[[4e5, -1.5e5], [4e5, 4.5e5], [6e5, 4.5e5], [6e5, -1.5e5], [4e5, -1.5e5]]] - @test b3.val == 0.0 - @test typeof(b3) == Subzero.CollisionBoundary{West, Float64} - @test b3.coords == [[[-2e5, -1.5e5], [-2e5, 4.5e5], [0.0, 4.5e5], [0.0, -1.5e5], [-2e5, -1.5e5]]] - @test b4.val == 0.0 - @test typeof(b4) == Subzero.PeriodicBoundary{South, Float64} - @test b4.coords == [[[-2e5, -1.5e5], [-2e5, 0.0], [6e5, 0.0], [6e5, -1.5e5], [-2e5, -1.5e5]]] - @test b5.u == 1.0 - @test b5.v == 2.0 - - # Test changing CollisionBoundary values and can't change other boundary types - b5.val = 1.0 - @test b5.val == 1.0 - @test_throws Exception b4.val = 1.0 - - # Creation of Float32 and Float64 Boundary - b32 = Subzero.OpenBoundary(Float32, North, g) - @test typeof(b32.val) == Float32 - @test typeof(b32) == Subzero.OpenBoundary{North, Float32} - b64 = Subzero.OpenBoundary(Float64, North, g) - @test typeof(b64.val) == Float64 - @test typeof(b64) == Subzero.OpenBoundary{North, Float64} - - # Periodic Compat - @test !Subzero.periodic_compat(b1, b2) - @test !Subzero.periodic_compat(b2, b1) - @test Subzero.periodic_compat(b1, b4) - @test Subzero.periodic_compat(b2, b3) - end - - @testset "Topography" begin - coords = [[[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]] - poly = Subzero.make_polygon(coords) - # Polygon Constructor - topo1 = Subzero.TopographyElement(poly) - @test topo1.coords == coords - @test topo1.centroid == [0.5, 0.5] - @test topo1.rmax == sqrt(0.5) - topo32 = Subzero.TopographyElement(Float32, poly) - @test typeof(topo32) == Subzero.TopographyElement{Float32} - @test typeof(topo32.coords) == Subzero.PolyVec{Float32} - # Coords Constructor - topo2 = Subzero.TopographyElement(Float64, coords) - @test topo2.coords == coords - @test topo2.centroid == [0.5, 0.5] - @test topo2.rmax == sqrt(0.5) - # Basic constructor - topo3 = TopographyElement( - Subzero.make_polygon(coords), - coords, - [0.5, 0.5], - sqrt(0.5), - ) - @test topo3.coords == coords - # check when radius is less than or equal to 0 - @test_throws ArgumentError TopographyElement( - Subzero.make_polygon(coords), - coords, - [0.5, 0.5], - -sqrt(0.5), - ) - - # Create field of topography - coords_w_hole = [ - [[0.5, 10.0], [0.5, 0.0], [10.0, 0.0], [10.0, 10.0], [0.5, 10.0]], - [[2.0, 8.0], [2.0, 4.0], [8.0, 4.0], [8.0, 8.0], [2.0, 8.0]] - ] - topo_field_64 = initialize_topography_field( - Float64, - [coords, coords_w_hole], - ) - @test length(topo_field_64) == 2 - @test typeof(topo_field_64) <: StructArray{TopographyElement{Float64}} - @test !Subzero.hashole(topo_field_64.coords[2]) - - topo_field_32 = initialize_topography_field( - Float32, - [coords, coords_w_hole], - ) - @test typeof(topo_field_32) <: StructArray{TopographyElement{Float32}} - - @test typeof(initialize_topography_field([coords, coords_w_hole])) <: - StructArray{TopographyElement{Float64}} - end - - @testset "Domain" begin - FT = Float64 - g = Subzero.RegRectilinearGrid( - (0, 4e5), - (0, 3e5), - 1e4, - 1e4, - ) - b1 = Subzero.PeriodicBoundary(North, g) - b2 = Subzero.OpenBoundary(East, g) - b3 = Subzero.CollisionBoundary(West, g) - b4 = Subzero.PeriodicBoundary(South, g) - topography = StructArray([Subzero.TopographyElement( - [[[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]], - )]) - # test basic domain with no topography - rdomain1 = Subzero.Domain(b1, b4, b2, b3) - @test rdomain1.north == b1 - @test rdomain1.south == b4 - @test rdomain1.east == b2 - @test rdomain1.west == b3 - @test isempty(rdomain1.topography) - # test basic domain with topography - rdomain2 = Subzero.Domain(b1, b4, b2, b3, topography) - @test isempty(rdomain1.topography) - # domain with wrong directions - @test_throws MethodError Subzero.Domain(b4, b2, b2, b3) - # domain with non-periodic - @test_throws ArgumentError Subzero.Domain( - b1, - Subzero.OpenBoundary(South, g), - b2, - b3, - ) - @test_throws ArgumentError Subzero.Domain( - b1, - b4, - b2, - Subzero.PeriodicBoundary(West, g), - ) - # domain with north < south - p_placeholder = GI.Polygon([[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)]]) - @test_throws ArgumentError Subzero.Domain( - b1, - Subzero.OpenBoundary( - South, - p_placeholder, - GI.coordinates(p_placeholder), - 6e5, - ), - b2, - b3, - ) - # domain with east < west - @test_throws ArgumentError Subzero.Domain( - b1, - b4, - b2, - Subzero.OpenBoundary( - West, - p_placeholder, - GI.coordinates(p_placeholder), - 6e5, - ), - ) - end - @testset "Model" begin # test domain in grid # test basic working model diff --git a/test/test_output.jl b/test/test_output.jl index 3ac285f..b5d4750 100644 --- a/test/test_output.jl +++ b/test/test_output.jl @@ -1,18 +1,13 @@ function test_basic_outputwriters() - grid = RegRectilinearGrid( - (-1e5, 1e5), - (-1e5, 1e5), - 1e4, - 1e4, - ) + grid = RegRectilinearGrid(; x0 = -1e5, xf = 1e5, y0 = -1e5, yf = 1e5, Δx = 1e4, Δy = 1e4) ocean = Ocean(grid, 0.0, 0.0, 0.0) atmos = Atmos(grid, 0.0, 0.0, 0.0) - domain = Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + domain = Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) floe_coords = [[[7.5e4, 7.5e4], [7.5e4, 9.5e4], [9.5e4, 9.5e4], [9.5e4, 7.5e4], [7.5e4, 7.5e4]]] diff --git a/test/test_physical_processes/test_collisions.jl b/test/test_physical_processes/test_collisions.jl index a2e7567..b9e9d92 100644 --- a/test/test_physical_processes/test_collisions.jl +++ b/test/test_physical_processes/test_collisions.jl @@ -9,33 +9,32 @@ # Grid setup Lx = 1e5 Ly = Lx - grid = RegRectilinearGrid((-Lx, Lx), (-Ly, Ly), 1e4, 1e4) + grid = RegRectilinearGrid(; x0 = -Lx, xf = Lx, y0 = -Ly, yf = Ly, Δx = 1e4, Δy = 1e4) # Boundary setup - p_n_bound = PeriodicBoundary(North, grid) - p_s_bound = PeriodicBoundary(South, grid) - p_e_bound = PeriodicBoundary(East, grid) - p_w_bound = PeriodicBoundary(West, grid) + p_n_bound = PeriodicBoundary(North; grid) + p_s_bound = PeriodicBoundary(South; grid) + p_e_bound = PeriodicBoundary(East; grid) + p_w_bound = PeriodicBoundary(West; grid) - c_n_bound = CollisionBoundary(North, grid) - c_s_bound = CollisionBoundary(South, grid) - c_e_bound = CollisionBoundary(East, grid) - c_w_bound = CollisionBoundary(West, grid) + c_n_bound = CollisionBoundary(North; grid) + c_s_bound = CollisionBoundary(South; grid) + c_e_bound = CollisionBoundary(East; grid) + c_w_bound = CollisionBoundary(West; grid) - o_n_bound = OpenBoundary(North, grid) - o_s_bound = OpenBoundary(South, grid) - o_e_bound = OpenBoundary(East, grid) - o_w_bound = OpenBoundary(West, grid) + o_n_bound = OpenBoundary(North; grid) + o_s_bound = OpenBoundary(South; grid) + o_e_bound = OpenBoundary(East; grid) + o_w_bound = OpenBoundary(West; grid) - topos = initialize_topography_field(FT, [[[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]]]) - # topos = StructArray([TopographyElement([[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]])]) + topos = initialize_topography_field(FT; coords = [[[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]]]) - topo_domain = Domain(p_n_bound, p_s_bound, c_e_bound, o_w_bound, topos) - collision_domain = Domain(c_n_bound, c_s_bound, c_e_bound, c_w_bound) - open_domain = Domain(o_n_bound, o_s_bound, o_e_bound, o_w_bound) - ew_periodic_domain = Domain(o_n_bound, o_s_bound, p_e_bound, p_w_bound) - ns_periodic_domain = Domain(p_n_bound, p_s_bound, o_e_bound, o_w_bound) - double_periodic_domain = Domain(p_n_bound, p_s_bound, p_e_bound, p_w_bound) + topo_domain = Domain(; north = p_n_bound, south = p_s_bound, east = c_e_bound, west = o_w_bound, topography = topos) + collision_domain = Domain(; north = c_n_bound, south = c_s_bound, east = c_e_bound, west = c_w_bound) + open_domain = Domain(; north = o_n_bound, south = o_s_bound, east = o_e_bound, west = o_w_bound) + ew_periodic_domain = Domain(; north = o_n_bound, south = o_s_bound, east = p_e_bound, west = p_w_bound) + ns_periodic_domain = Domain(; north = p_n_bound, south = p_s_bound, east = o_e_bound, west = o_w_bound) + double_periodic_domain = Domain(; north = p_n_bound, south = p_s_bound, east = p_e_bound, west = p_w_bound) @testset "Floe-Floe Interactions" begin hmean = 0.25 @@ -186,30 +185,6 @@ Subzero.floe_domain_interaction!(corner_floe, collision_domain, consts, Δt, max_overlap) @test all(corner_floe.interactions[:, xforce] .<= 0) @test all(corner_floe.interactions[:, yforce] .<= 0) - - # Test compression boundaries movement - TODO: Move these to different file - nc_boundary = MovingBoundary(North, grid, 0.0, -0.1) - nc_coords = deepcopy(nc_boundary.coords) - sc_boundary = MovingBoundary(South, grid, 0.0, 0.1) - sc_coords = deepcopy(sc_boundary.coords) - ec_boundary = MovingBoundary(East, grid, 0.1, 0.0) - ec_coords = deepcopy(ec_boundary.coords) - wc_boundary = MovingBoundary(West, grid, 0.1, 0.0) - wc_coords = deepcopy(wc_boundary.coords) - cdomain = Domain(nc_boundary, sc_boundary, ec_boundary, wc_boundary) - Subzero.update_boundaries!(cdomain, 10) - Subzero.translate!(nc_coords, 0, -1) - @test nc_coords == nc_boundary.coords - @test nc_boundary.val == 1e5 - 1 - Subzero.translate!(sc_coords, 0, 1) - @test sc_coords == sc_boundary.coords - @test sc_boundary.val == -1e5 + 1 - Subzero.translate!(ec_coords, 1, 0) - @test ec_coords == ec_boundary.coords - @test ec_boundary.val == 1e5 + 1 - Subzero.translate!(wc_coords, 1, 0) - @test wc_coords == wc_boundary.coords - @test wc_boundary.val == -1e5 + 1 end @testset "Add Ghosts" begin diff --git a/test/test_physical_processes/test_coupling.jl b/test/test_physical_processes/test_coupling.jl index e0ee7c8..ee48b20 100644 --- a/test/test_physical_processes/test_coupling.jl +++ b/test/test_physical_processes/test_coupling.jl @@ -163,12 +163,7 @@ end @testset "Coupling Helper Functions" begin - grid = Subzero.RegRectilinearGrid( - (-10, 10), - (-8, 8), - 2, - 4, - ) + grid = Subzero.RegRectilinearGrid(; x0 = -10, xf = 10, y0 = -8, yf = 8, Δx = 2, Δy = 4) # Test find_cell_indices xpoints = [-10.5, -10, -10, -6.5, -6, -4, 10, 10.5, 12] ypoints = [0.0, 6.0, -8.0, 4.5, 0.0, 5.0, -8.0, 0.0, 0.0] @@ -183,8 +178,8 @@ end # Test filter_oob_points - open_bound = Subzero.OpenBoundary(East, grid) - periodic_bound = Subzero.PeriodicBoundary(East, grid) + open_bound = Subzero.OpenBoundary(East; grid) + periodic_bound = Subzero.PeriodicBoundary(East; grid) x = [-12, -10, -8, -6, 0, 4, 4, 10, 12, 12] y = [5, -6, 4, 10, -10, 8, -8, -6, 4, 10] open_open_answers = [false, true, true, false, false, true, true, true, false, false] @@ -469,19 +464,14 @@ @testset "OA Forcings" begin # # set up model and floe FT = Float64 - grid = Subzero.RegRectilinearGrid( - (-1e5, 1e5), - (-1e5, 1e5), - 1e4, - 1e4, - ) + grid = Subzero.RegRectilinearGrid(; x0 = -1e5, xf = 1e5, y0 = -1e5, yf = 1e5, Δx = 1e4, Δy = 1e4) zonal_ocean = Subzero.Ocean(grid, 1.0, 0.0, 0.0) zero_atmos = Subzero.Atmos(grid, 0.0, 0.0, -20.0) - domain = Subzero.Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), + domain = Subzero.Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = CollisionBoundary(East; grid), + west = CollisionBoundary(West; grid), ) floe = Subzero.Floe( [[ diff --git a/test/test_physical_processes/test_ridge_raft.jl b/test/test_physical_processes/test_ridge_raft.jl index aa6772d..bfedeb9 100644 --- a/test/test_physical_processes/test_ridge_raft.jl +++ b/test/test_physical_processes/test_ridge_raft.jl @@ -254,29 +254,24 @@ end # Setup for tests - grid = RegRectilinearGrid( - (0, 1e5), - (0, 1e5), - 1e4, - 1e4, - ) + grid = RegRectilinearGrid(; x0 = 0, xf = 1e5, y0 = 0, yf = 1e5, Δx = 1e4, Δy = 1e4) topo_coords = [[[5e4, 5e4], [5e4, 7e4], [7e4, 7e4], [7e4, 5e4], [5e4, 5e4]]] - collision_domain = Subzero.Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - CollisionBoundary(East, grid), - CollisionBoundary(West, grid), - initialize_topography_field([topo_coords]) + collision_domain = Subzero.Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = CollisionBoundary(East; grid), + west = CollisionBoundary(West; grid), + topography = initialize_topography_field(; coords = [topo_coords]) ) - periodic_domain = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), + periodic_domain = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), ) boundary_poly = GO.UnionIntersectingPolygons()(Subzero.make_multipolygon([ - collision_domain.north.coords, collision_domain.south.coords, - collision_domain.east.coords, collision_domain.west.coords])) + collision_domain.north.poly, collision_domain.south.poly, + collision_domain.east.poly, collision_domain.west.poly])) topo_poly = Subzero.make_polygon(topo_coords) consts = Constants() diff --git a/test/test_physical_processes/test_simplification.jl b/test/test_physical_processes/test_simplification.jl index 5f4af56..0686a59 100644 --- a/test/test_physical_processes/test_simplification.jl +++ b/test/test_physical_processes/test_simplification.jl @@ -2,17 +2,12 @@ FT = Float64 Δt = 10 @testset "Dissolve Floes" begin - grid = RegRectilinearGrid( - (-1e5, 1e5), - (0.0, 1e5), - 1e4, - 1e4, - ) - domain = Subzero.Domain( - CollisionBoundary(North, grid), - CollisionBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), + grid = RegRectilinearGrid(; x0 = -1e5, xf = 1e5, y0 = 0.0, yf = 1e5, Δx = 1e4, Δy = 1e4) + domain = Subzero.Domain(; + north = CollisionBoundary(North; grid), + south = CollisionBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), ) height = 0.25 ρi = 920.0 @@ -210,17 +205,12 @@ @test f3.status.tag == Subzero.active # Test overall fuse floe functionality with set of 4 floes - grid = RegRectilinearGrid( - (-2.5e4, 1e5), - (-2.5e4, 1e5), - 1e4, - 1e4, - ) - open_domain_no_topo = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + grid = RegRectilinearGrid(; x0 = -2.5e4, xf = 1e5, y0 = -2.5e4, yf = 1e5, Δx = 1e4, Δy = 1e4) + open_domain_no_topo = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) coords1 = [[ # large floe [0.0, 0.0], @@ -287,17 +277,12 @@ @test floe_arr.area[3] == floe3_area # small floes fused into floe 1 end @testset "Smooth Floes" begin - grid = RegRectilinearGrid( - (-2.5e4, 1e5), - (-2.5e4, 1e5), - 1e4, - 1e4, - ) - open_domain_no_topo = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + grid = RegRectilinearGrid(; x0 = -2.5e4, xf = 1e5, y0 = -2.5e4, yf = 1e5, Δx = 1e4, Δy = 1e4) + open_domain_no_topo = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) # Create complex floes file = jldopen("inputs/floe_shapes.jld2", "r") @@ -429,22 +414,13 @@ # Two floes overlap, and one is cut into two pieces by topography floe_set2 = floe_arr[1:2] - open_domain_with_topo = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), - StructVector( - [TopographyElement( - [[ - [0.0, 1.05e4], - [0.0, 1.15e4], - [3e3, 1.15e4], - [3e3, 1.05e4], - [0.0, 1.05e4], - ]], - )], - ), + open_domain_with_topo = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), + topography = initialize_topography_field(; coords = [[[ + [0.0, 1.05e4], [0.0, 1.15e4], [3e3, 1.15e4], [3e3, 1.05e4], [0.0, 1.05e4]]]]) ) og_f1_area = floe_set2.area[1] total_mass = sum(floe_set2.mass) @@ -474,17 +450,13 @@ @test floe_set2.status[2].fuse_idx == [1] end @testset "Remove Floes" begin - grid = RegRectilinearGrid( - (-2.5e4, 1e5), - (-2.5e4, 1e5), - 1e4, - 1e4, - ) - open_domain_no_topo = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + grid = RegRectilinearGrid(; x0 = -2.5e4, xf = 1e5, y0 = -2.5e4, yf = 1e5, Δx = 1e4, Δy = 1e4) + + open_domain_no_topo = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) coords1 = [[ # large floe [0.0, 0.0], diff --git a/test/test_physical_processes/test_welding.jl b/test/test_physical_processes/test_welding.jl index 000d695..908eac8 100644 --- a/test/test_physical_processes/test_welding.jl +++ b/test/test_physical_processes/test_welding.jl @@ -1,29 +1,24 @@ @testset "Bin floes" begin Δt = 10 - grid = RegRectilinearGrid( - (0, 1e5), - (0, 1e5), - 1e4, - 1e4, - ) + grid = RegRectilinearGrid(; x0 = 0.0, xf = 1e5, y0 = 0.0, yf = 1e5, Δx = 1e4, Δy = 1e4) - open_domain = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + open_domain = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) - periodic_domain = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - PeriodicBoundary(East, grid), - PeriodicBoundary(West, grid), + periodic_domain = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = PeriodicBoundary(East; grid), + west = PeriodicBoundary(West; grid), ) - half_open_periodic_domain = Subzero.Domain( - PeriodicBoundary(North, grid), - PeriodicBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + half_open_periodic_domain = Subzero.Domain(; + north = PeriodicBoundary(North; grid), + south = PeriodicBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) coords = [ @@ -134,17 +129,13 @@ end @testset "Weld floes" begin Δt = 10 - grid = RegRectilinearGrid( - (0, 1e5), - (0, 1e5), - 1e4, - 1e4, - ) - periodic_domain = Subzero.Domain( - OpenBoundary(North, grid), - OpenBoundary(South, grid), - OpenBoundary(East, grid), - OpenBoundary(West, grid), + grid = RegRectilinearGrid(; x0 = 0.0, xf = 1e5, y0 = 0.0, yf = 1e5, Δx = 1e4, Δy = 1e4) + + periodic_domain = Subzero.Domain(; + north = OpenBoundary(North; grid), + south = OpenBoundary(South; grid), + east = OpenBoundary(East; grid), + west = OpenBoundary(West; grid), ) consts = Constants() floe_settings = FloeSettings() diff --git a/test/test_simulation_components/domain_components/boundaries.jl b/test/test_simulation_components/domain_components/boundaries.jl new file mode 100644 index 0000000..496b6fc --- /dev/null +++ b/test/test_simulation_components/domain_components/boundaries.jl @@ -0,0 +1,132 @@ +using Test, Subzero +import GeometryOps as GO +import GeometryOps.GeoInterface as GI + +@testset "Directions" begin + FT = Float64 + x0, xf, y0, yf = 0.0, 1e5, -5e4, 5e4 + Δx, Δy = (xf - x0) / 2, (yf - y0) / 2 + # North boundary + n_poly, n_val = Subzero._boundary_info_from_extent(North, FT, x0, xf, y0, yf) + n_point_set = Set(GI.getpoint(n_poly)) + @test issetequal(n_point_set, Set(((x0 - Δx, yf), (x0 - Δx, yf + Δy), (xf + Δx, yf + Δy), (xf + Δx, yf)))) + @test n_val == yf + # South boundary + s_poly, s_val = Subzero._boundary_info_from_extent(South, FT, x0, xf, y0, yf) + s_point_set = Set(GI.getpoint(s_poly)) + @test issetequal(s_point_set, Set(((x0 - Δx, y0), (x0 - Δx, y0 - Δy), (xf + Δx, y0 - Δy), (xf + Δx, y0)))) + @test s_val == y0 + # East boundary + e_poly, e_val = Subzero._boundary_info_from_extent(East, FT, x0, xf, y0, yf) + e_point_set = Set(GI.getpoint(e_poly)) + @test issetequal(e_point_set, Set(((xf, y0 - Δy), (xf, yf + Δy), (xf + Δx, y0 - Δy), (xf + Δx, yf + Δy)))) + @test e_val == xf + # West boundary + w_poly, w_val = Subzero._boundary_info_from_extent(West, FT, x0, xf, y0, yf) + w_point_set = Set(GI.getpoint(w_poly)) + @test issetequal(w_point_set, Set(((x0, y0 - Δy), (x0, yf + Δy), (x0 - Δx, y0 - Δy), (x0 - Δx, yf + Δy)))) + @test w_val == x0 +end + +@testset "Boundaries" begin + x0, xf, y0, yf = 0, 4e5, 0, 3e5 + grid = RegRectilinearGrid(; x0, xf, y0, yf, Δx = 1e4, Δy = 1e4) + # Check Northern Boundary Values + n_p = PeriodicBoundary(North; grid) + n_o = OpenBoundary(North; x0, xf, y0, yf) + # ensure val and poly fields are the same (and correct) regardless of constructor branch + @test n_p.val == n_o.val == 3e5 + @test GO.equals(n_p.poly, n_o.poly) + @test GO.area(n_p.poly) == 8e5 * 1.5e5 + north_points = Set(((-2e5, 3e5), (-2e5, 4.5e5), (6e5, 4.5e5), (6e5, 3e5))) + @test issetequal(Set(GI.getpoint(n_p.poly)), north_points) + # check type correctness + @test typeof(n_p) == PeriodicBoundary{North, Float64} + @test typeof(n_o) == OpenBoundary{North, Float64} + + # Check Eastern Boundary Values + e_o = OpenBoundary(East; grid) + e_c = CollisionBoundary(East; x0, xf, y0, yf) + # ensure val and poly fields are the same (and correct) regardless of constructor branch + @test e_o.val == e_c.val == 4e5 + @test GO.equals(e_o.poly, e_c.poly) + @test GO.area(e_o.poly) == 2e5 * 6e5 + east_points = Set(((4e5, -1.5e5), (4e5, 4.5e5), (6e5, 4.5e5), (6e5, -1.5e5))) + @test issetequal(Set(GI.getpoint(e_o.poly)), east_points) + # check type correctness + @test typeof(e_o) == OpenBoundary{East, Float64} + @test typeof(e_c) == CollisionBoundary{East, Float64} + + # Check Southern Boundary Values + s_p = PeriodicBoundary(South, Float64; grid) + s_m = MovingBoundary(South, Float64; grid, u = 1.0, v = 2.0) + # ensure val and poly fields are the same (and correct) regardless of constructor branch + @test s_p.val == s_m.val == 0 + @test GO.equals(s_p.poly, s_m.poly) + @test GO.area(s_p.poly) == 8e5 * 1.5e5 + south_points = Set(((-2e5, -1.5e5), (-2e5, 0.0), (6e5, 0.0), (6e5, -1.5e5))) + s_p_points = Set(GI.getpoint(s_p.poly)) + @test issetequal(s_p_points, south_points) + # check type correctness + @test typeof(s_p) == PeriodicBoundary{South, Float64} + @test typeof(s_m) == MovingBoundary{South, Float64} + @test eltype(first(s_p_points)) == typeof(s_p.val) == Float64 + + # Check Western Boundary Values + w_c = CollisionBoundary(West, Float32; grid) + w_m = MovingBoundary(West, Float32; x0, xf, y0, yf, u = 0.1, v = -0.1) + # ensure val and poly fields are the same (and correct) regardless of constructor branch + @test w_c.val == w_m.val == 0 + @test GO.equals(w_c.poly, w_m.poly) + @test GO.area(w_c.poly, Float32) ≈ 2e5 * 6e5 + west_points = Set(((-2e5, -1.5e5), (-2e5, 4.5e5), (0.0, 4.5e5), (0.0, -1.5e5))) + w_c_points = Set(GI.getpoint(w_c.poly)) + @test issetequal(w_c_points, west_points) + # check type correctness + @test typeof(w_c) == CollisionBoundary{West, Float32} + @test typeof(w_m) == MovingBoundary{West, Float32} + @test eltype(first(w_c_points)) == typeof(w_c.val) == Float32 + + # Check _get_velocity boundary dispatch + x, y = 0.0, 0.0 # (x, y) point does not matter as boundaries don't rotate + @test Subzero._get_velocity(s_m, x, y) == (s_m.u, s_m.v) == (1.0, 2.0) + @test Subzero._get_velocity(w_m, x, y) == (w_m.u, w_m.v) == (0.1f0, -0.1f0) + @test Subzero._get_velocity(n_o, x, y) == (0.0, 0.0) + @test Subzero._get_velocity(e_c, x, y) == (0.0, 0.0) + @test Subzero._get_velocity(s_p, x, y) == (0.0, 0.0) + + # Check _update_boundary! boundary dispatch + FT, Δt = Float64, 20 + # check that OpenBoundary does not move + curr_val, curr_poly = n_o.val, n_o.poly + Subzero._update_boundary!(n_o, Δt) + @test n_o.val == curr_val && GO.equals(n_o.poly, curr_poly) + # check that CollisionBoundary does not move + curr_val, curr_poly = e_c.val, e_c.poly + Subzero._update_boundary!(e_c, Δt) + @test e_c.val == curr_val && GO.equals(e_c.poly, curr_poly) + # check that PeriodicBoundary does not move + curr_val, curr_poly = s_p.val, s_p.poly + Subzero._update_boundary!(s_p, Δt) + @test s_p.val == curr_val && GO.equals(s_p.poly, curr_poly) + # check that Western MovingBoundary DOES move + curr_val, curr_poly = w_m.val, w_m.poly + Subzero._update_boundary!(w_m, Δt) + Δx = Δt * w_m.u + @test w_m.u != 0 + @test w_m.val == curr_val + Δx + @test GO.equals(w_m.poly, Subzero._translate_poly(FT, curr_poly, Δx, 0.0)) + # check that Southern MovingBoundary DOES move + curr_val, curr_poly = s_m.val, s_m.poly + Subzero._update_boundary!(s_m, Δt) + Δy = Δt * s_m.v + @test s_m.v != 0 + @test s_m.val == curr_val + Δy + @test GO.equals(s_m.poly, Subzero._translate_poly(FT, curr_poly, 0.0, Δy)) + + # Check Periodic Compatability + @test !Subzero._periodic_compat(n_p, s_m) + @test !Subzero._periodic_compat(n_o, s_p) + @test Subzero._periodic_compat(n_p, s_p) + @test Subzero._periodic_compat(e_c, w_m) +end diff --git a/test/test_simulation_components/domain_components/domains.jl b/test/test_simulation_components/domain_components/domains.jl new file mode 100644 index 0000000..8b3cd10 --- /dev/null +++ b/test/test_simulation_components/domain_components/domains.jl @@ -0,0 +1,65 @@ +using Test, Subzero +import GeometryOps as GO +import GeometryOps.GeoInterface as GI + +@testset "Domain" begin + grid = RegRectilinearGrid(; x0 = 0, xf = 4e5, y0 = 0, yf = 3e5, Δx = 1e4, Δy = 1e4) + north_periodic = PeriodicBoundary(North; grid) + north_collision = CollisionBoundary(North; grid) + south_periodic = PeriodicBoundary(South; grid) + south_collision = CollisionBoundary(South; grid) + east_periodic = PeriodicBoundary(East; grid) + east_collision = CollisionBoundary(East; grid) + west_periodic = PeriodicBoundary(West; grid) + west_collision = CollisionBoundary(West; grid) + + # test basic periodic domain with no topography (creates default topography) + domain1 = Domain(; north = north_periodic, south = south_periodic, east = east_periodic, west = west_periodic) + @test typeof(domain1) <: Domain{Float64, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary} + @test domain1.north == north_periodic + @test domain1.south == south_periodic + @test domain1.east == east_periodic + @test domain1.west == west_periodic + @test typeof(domain1.topography) <: Subzero.TopographyField{Float64} + @test isempty(domain1.topography) + + # test basic non-periodic domain with topography + topography = initialize_topography_field(; coords = [[[(1e4, 1e4), (1e4, 2e4), (2e4, 2e4), (2e4, 1e4), (1e4, 1e4)]]]) + domain2 = Domain(; topography, north = north_collision, south = south_collision, east = east_collision, west = west_collision) + @test domain2.north == north_collision + @test domain2.south == south_collision + @test domain2.east == east_collision + @test domain2.west == west_collision + @test typeof(domain2.topography) <: Subzero.TopographyField{Float64} + @test length(domain2.topography) == 1 + @test GO.area(domain2.topography[1].poly) == 1e8 + + # test error with wrong directions paired with wrong keywords --> error + @test_throws MethodError Domain(; south = north_collision, north = south_collision, east = east_collision, west = west_collision) + @test_throws MethodError Domain(; north = north_collision, south = south_collision, west = east_collision, east = west_collision) + + # domain with non-compatible periodic pairs --> error + @test_throws ArgumentError Domain(; north = north_collision, south = south_periodic, east = east_collision, west = west_collision) + @test_throws ArgumentError Domain(; north = north_collision, south = south_collision, east = east_periodic, west = west_collision) + + # domain with north.val < south.val --> error + north_low = PeriodicBoundary(North; x0 = 0.0, xf = 4e5, y0 = -4e5, yf = 0.0) + south_high = PeriodicBoundary(South; x0 = 0.0, xf = 4e5, y0 = 4e5, yf = 8e5) + @test_throws ArgumentError Domain(; north = north_low, south = south_high, east = east_collision, west = west_collision) + + # domain with east < west --> error + east_low = PeriodicBoundary(East; x0 = -4e5, xf = 0.0, y0 = 0.0, yf = 4e5) + west_high = PeriodicBoundary(West; x0 = 4e5, xf = 8e5, y0 = 0.0, yf = 4e5) + @test_throws ArgumentError Domain(; north = north_collision, south = south_collision, east = east_low, west = west_high) + + # test Float32 Domain + north_periodic_32 = PeriodicBoundary(North, Float32; grid) + south_periodic_32 = PeriodicBoundary(South, Float32; grid) + east_periodic_32 = PeriodicBoundary(East, Float32; grid) + west_periodic_32 = PeriodicBoundary(West, Float32; grid) + domain3 = Domain(; north = north_periodic_32, south = south_periodic_32, east = east_periodic_32, west = west_periodic_32) + @test typeof(domain3) <: Domain{Float32, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary} + + # domain with mixed types --> error + @test_throws MethodError Domain(; north = north_periodic, south = south_periodic_32, east = east_periodic_32, west = west_periodic_32) +end \ No newline at end of file diff --git a/test/test_simulation_components/domain_components/topography.jl b/test/test_simulation_components/domain_components/topography.jl new file mode 100644 index 0000000..5462e56 --- /dev/null +++ b/test/test_simulation_components/domain_components/topography.jl @@ -0,0 +1,93 @@ +using Test, Subzero +import GeometryOps as GO +import GeometryOps.GeoInterface as GI +import LibGEOS as LG + +@testset "Topography" begin + # Test basic polygon + poly1 = GI.Polygon([[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]]) + topo1_Float64 = Subzero.TopographyElement(; poly = poly1) + # check types + @test typeof(topo1_Float64) == TopographyElement{Float64} + @test GI.x(topo1_Float64.centroid) isa Float64 && GI.y(topo1_Float64.centroid) isa Float64 + @test topo1_Float64.rmax isa Float64 + # check values + @test GO.equals(poly1, topo1_Float64.poly) + @test all(topo1_Float64.centroid .≈ (0.5, 0.5)) + @test topo1_Float64.rmax ≈ sqrt(0.5^2 + 0.5^2) + + # Test basic polygon with Float32 + topo1_Float32 = Subzero.TopographyElement(Float32; poly = poly1) + # check types + @test typeof(topo1_Float32) == TopographyElement{Float32} + @test GI.x(topo1_Float32.centroid) isa Float32 && GI.y(topo1_Float32.centroid) isa Float32 + @test topo1_Float32.rmax isa Float32 + + # Test non-closed ring polygon + poly2 = GI.Polygon([[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)]]) + topo2 = Subzero.TopographyElement(; poly = poly2) + # check close polygon + @test GI.npoint(topo2.poly) == 5 + topo2_ring = GI.getexterior(topo2.poly) + @test (GI.getpoint(topo2_ring, 1) == GI.getpoint(topo2_ring, 5)) + + # Test polygon with hole + poly3 = GI.Polygon([ + [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)], + [(0.2, 0.2), (0.3, 0.2), (0.3, 0.3), (0.2, 0.3), (0.2, 0.2)], + ]) + topo3 = Subzero.TopographyElement(; poly = poly3) + # check polygon has no holes + @test GI.nhole(topo3.poly) == 0 + + # Test TopographyElement _get_velocity + x, y = 0.5, 0.5 + @test Subzero._get_velocity(topo1_Float64, x, y) == (0.0, 0.0) + + # Test TopographyElement _normal_direction_correct! + forces, points = ones(2, 2), ones(2, 2) + Subzero._normal_direction_correct!(forces, points, topo1_Float64) + @test forces == ones(2, 2) +end + +@testset "TopographyField" begin + coords1 = [[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]] + coords2 = [[[0.5, 0.0], [0.5, 1.0], [1.5, 1.0], [1.5, 0.0], [0.5, 0.0]]] + poly1 = GI.Polygon(coords1) + poly2 = LG.Polygon(coords2) + + # Topography field of 2 identical polygons is reduced to one element + topo_field_coords_64 = initialize_topography_field(; coords = [coords1, coords1]) + topo_field_polys_64 = initialize_topography_field(; polys = [poly1, poly1]) + # test field types + @test typeof(topo_field_coords_64) == typeof(topo_field_polys_64) + @test typeof(topo_field_coords_64) <: Subzero.TopographyField{Float64} + # test polygons are the same and equal to poly1 + @test length(topo_field_coords_64) == length(topo_field_polys_64) == 1 + @test GO.equals(topo_field_coords_64.poly[1], poly1) + @test GO.equals(topo_field_polys_64.poly[1], poly1) + # test centroid and rmax + @test topo_field_coords_64.centroid[1] == topo_field_polys_64.centroid[1] == [0.5, 0.5] + @test topo_field_coords_64.rmax[1] == topo_field_polys_64.rmax[1] == sqrt(0.5^2 + 0.5^2) + + # Topography field of 2 overlapping polygons is make into two non-intersecting elements + topo_field_coords_32 = initialize_topography_field(Float32; coords = [coords1, coords2]) + topo_field_polys_32 = initialize_topography_field(Float32; polys = [poly1, poly2]) + # test field types + @test typeof(topo_field_coords_32) == typeof(topo_field_polys_32) + @test typeof(topo_field_polys_32) <: Subzero.TopographyField{Float32} + # test two constructor methods produce the same two polygons + @test length(topo_field_coords_32) == length(topo_field_polys_32) == 2 + @test GO.equals(topo_field_coords_32.poly[1], topo_field_polys_32.poly[1]) + @test GO.equals(topo_field_coords_32.poly[2], topo_field_polys_32.poly[2]) + # test that the two polygons do not overlap + (p1_x0, p1_xf), (p1_y0, p1_yf) = GI.extent(topo_field_coords_32.poly[1]) + @test (p1_x0, p1_xf) == (0.0, 0.5) && (p1_y0, p1_yf) == (0.0, 1.0) + (p2_x0, p2_xf), (p2_y0, p2_yf) = GI.extent(topo_field_coords_32.poly[2]) + @test (p2_x0, p2_xf) == (0.5, 1.5) && (p2_y0, p2_yf) == (0.0, 1.0) + # test centroid and rmax + @test topo_field_coords_32.centroid[1] == topo_field_coords_32.centroid[1] == [0.25, 0.5] + @test topo_field_coords_32.centroid[2] == topo_field_coords_32.centroid[2] == [1.0, 0.5] + @test topo_field_coords_32.rmax[1] ≈ topo_field_coords_32.rmax[1] ≈ sqrt(0.25^2 + 0.5^2) + @test topo_field_coords_32.rmax[2] ≈ topo_field_coords_32.rmax[2] ≈ sqrt(0.5^2 + 0.5^2) +end diff --git a/test/test_simulation_components/grids.jl b/test/test_simulation_components/grids.jl new file mode 100644 index 0000000..8401ba9 --- /dev/null +++ b/test/test_simulation_components/grids.jl @@ -0,0 +1,50 @@ +using Test, Subzero + +@testset "CellFloes" begin + # Empty CellFloes + c1 = Subzero.CellFloes() + @test isempty(c1.floeidx) && isempty(c1.Δx) && isempty(c1.Δy) + @test typeof(c1) == Subzero.CellFloes{Float64} + # Non-empty CellFloes + c2 = Subzero.CellFloes(Float32; floeidx = [1, 2], Δx = [0.0, 0.0], Δy = [4.0, 4.0]) + @test length(c2.floeidx) == length(c2.Δx) == length(c2.Δy) == 2 + @test typeof(c2) == Subzero.CellFloes{Float32} + # CellFloes with incorrect dimensions + @test_throws AssertionError Subzero.CellFloes(Float32; floeidx = [1], Δx = [0.0, 0.0], Δy = [4.0, 4.0]) + # Empty CellFloes + empty!(c2) + @test isempty(c2.floeidx) && isempty(c2.Δx) && isempty(c2.Δy) +end + +@testset "RegRectilinearGrid" begin + # Non-square grid using constructor with Δx and Δy + g1 = Subzero.RegRectilinearGrid(Float64; x0 = -10, xf = 10, y0 = -8, yf = 8, Δx = 2, Δy = 4) + @test g1.Nx == 10 + @test g1.Ny == 4 + @test g1.x0 == -10 + @test g1.xf == 10 + @test g1.y0 == -8 + @test g1.yf == 8 + @test size(g1.floe_locations) == (11, 5) + @test typeof(g1) == Subzero.RegRectilinearGrid{Float64} + + # Non-square grid using constructor with Nx and Ny + g2 = Subzero.RegRectilinearGrid(Float32; x0 = -10, xf = 10, y0 = -8, yf = 8, Nx = 10, Ny = 4) + @test g2.x0 == -10 + @test g2.xf == 10 + @test g2.y0 == -8 + @test g2.yf == 8 + @test g2.Δx == 2 + @test g2.Δy == 4 + @test size(g2.floe_locations) == (11, 5) + @test typeof(g2) == Subzero.RegRectilinearGrid{Float32} + + # Test grid clipping to fit requested Δx and Δy when they divide evenly into (xf - x0) and (yf - y0) + g3 = Subzero.RegRectilinearGrid(; x0 = 0, xf = 20, y0 = 0, yf = 16, Δx = 3, Δy = 3) + @test g3.xf == 18 + @test g3.yf == 15 + @test g3.Nx == 6 + @test g3.Ny == 5 + @test size(g3.floe_locations) == (7, 6) + @test typeof(g3) == Subzero.RegRectilinearGrid{Float64} +end \ No newline at end of file