From 9a195cdb690090e5bbb11c0fe9a5c5ddb20f8fec Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Tue, 23 Apr 2024 11:21:08 -0700 Subject: [PATCH] Add floe creation bounds (#89) * New example file * Add floe creation bounds within domain * Add comment --- examples/forcing_contained_floes.jl | 106 ++++++++++++++++++++++++++++ src/simulation_components/floe.jl | 38 +++++----- 2 files changed, 126 insertions(+), 18 deletions(-) create mode 100644 examples/forcing_contained_floes.jl diff --git a/examples/forcing_contained_floes.jl b/examples/forcing_contained_floes.jl new file mode 100644 index 0000000..df2d169 --- /dev/null +++ b/examples/forcing_contained_floes.jl @@ -0,0 +1,106 @@ +using JLD2, Random, Statistics, Subzero + +# User Inputs +const FT = Float64 +const Lx = 1e5 +const Ly = 1e5 +const Δgrid = 2e3 +const hmean = 0.25 +const Δh = 0.0 +const Δt = 20 + +# Model instantiation +grid = RegRectilinearGrid( + (0.0, Lx), + (0.0, Ly), + Δgrid, + Δgrid, +) +# Set ocean u velocities +ocean_uvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1)) +for i in CartesianIndices(ocean_uvels) + r, c = Tuple(i) + if r ≤ 5 + ocean_uvels[i] = 0.2 + elseif r ≥ grid.Nx - 4 + ocean_uvels[i] = -0.2 + end +end +ocean_uvels[20:40, 20:30] .= 0.15 +# Set ocean v velocities +ocean_vvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1)) +for i in CartesianIndices(ocean_vvels) + r, c = Tuple(i) + if c ≤ 5 + ocean_vvels[i] = 0.2 + elseif c ≥ grid.Ny - 4 + ocean_vvels[i] = -0.2 + end +end +# Create ocean +ocean = Ocean( + ocean_uvels, + ocean_vvels, + zeros(FT, (grid.Nx + 1, grid.Ny + 1)), +) + +# Create atmosphere +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) + +domain = Domain(nboundary, sboundary, eboundary, wboundary) + +floe_settings = FloeSettings( + subfloe_point_generator = SubGridPointsGenerator(grid, 2), +) +# Floe creation - bound floes within smaller part of the domain +floe_bounds = [[[0.1Lx, 0.1Ly], [0.1Lx, 0.9Ly], [0.9Lx, 0.9Ly], [0.9Lx, 0.1Ly], [0.1Lx, 0.1Ly]]] +floe_arr = initialize_floe_field( + FT, + 300, + [0.4], + domain, + hmean, + Δh; + floe_bounds = floe_bounds, + rng = Xoshiro(1), + floe_settings = floe_settings +) + +# Model creation +model = Model(grid, ocean, atmos, domain, floe_arr) + +# Simulation setup +modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area))) +consts = Constants(E = modulus) + +# Run simulation +run_time!(simulation) = @time run!(simulation) +dir = "output/contained" + +# Output setup +initwriter = InitialStateOutputWriter(dir = dir, overwrite = true) +floewriter = FloeOutputWriter(50, dir = dir, overwrite = true) +writers = OutputWriters(initwriter, floewriter) +simulation = Simulation( + model = model, + consts = consts, + Δt = Δt, + nΔt = 15000, + verbose = true, + writers = writers, + rng = Xoshiro(1), +) +run_time!(simulation) + +plot_sim( + joinpath(dir, "floes.jld2"), + joinpath(dir, "initial_state.jld2"), + Δt, + joinpath(dir, "contained.mp4"), +) \ No newline at end of file diff --git a/src/simulation_components/floe.jl b/src/simulation_components/floe.jl index ab0908c..9c944f8 100644 --- a/src/simulation_components/floe.jl +++ b/src/simulation_components/floe.jl @@ -606,9 +606,12 @@ Inputs: hmean average floe height Δh height range - floes will range in height from hmean - Δh to hmean + Δh - floe_settings settings needed to initialize floes - rng random number generator to generate random floe - attributes - default uses Xoshiro256++ + floe_bounds coordinates of boundary within which to populate floes. This + can be smaller that the domain, but will be limited to open space + within the domain + floe_settings settings needed to initialize floes + rng random number generator to generate random floe + attributes - default uses Xoshiro256++ Output: floe_arr list of floes created using Voronoi Tesselation of the domain with given concentrations. @@ -620,31 +623,30 @@ function initialize_floe_field( domain, hmean, Δh; + floe_bounds = rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val), floe_settings = FloeSettings(min_floe_area = 0.0), rng = Xoshiro(), ) where {FT <: AbstractFloat} floe_arr = StructArray{Floe{FT}}(undef, 0) - # Split domain into cells with given concentrations - nrows, ncols = size(concentrations[:, :]) - Lx = domain.east.val - domain.west.val - Ly = domain.north.val - domain.south.val - rowlen = Ly / nrows - collen = Lx / ncols - # Availible space in whole domain - open_water = LG.Polygon(rect_coords( - domain.west.val, - domain.east.val, - domain.south.val, - domain.north.val - )) + # Availible space in domain + open_water = LG.Polygon(floe_bounds) + open_water = LG.intersection(open_water, LG.Polygon(rect_coords(domain.west.val, domain.east.val, domain.south.val, domain.north.val))) if !isempty(domain.topography) open_water = LG.difference( open_water, LG.MultiPolygon(domain.topography.coords) ) end + (bounds_xmin, bounds_xmax), (bounds_ymin, bounds_ymax) = GeometryBasics.GeoInterface.extent(open_water) open_water_area = LG.area(open_water) + # Split domain into cells with given concentrations + nrows, ncols = size(concentrations[:, :]) + Lx = bounds_xmax - bounds_xmin + Ly = bounds_ymax - bounds_ymin + rowlen = Ly / nrows + collen = Lx / ncols + # Loop over cells for j in range(1, ncols) for i in range(1, nrows) @@ -652,8 +654,8 @@ function initialize_floe_field( if c > 0 c = c > 1 ? 1 : c # Grid cell bounds - xmin = domain.west.val + collen * (j - 1) - ymin = domain.south.val + rowlen * (i - 1) + xmin = bounds_xmin + collen * (j - 1) # TODO: change this + ymin = bounds_ymin + rowlen * (i - 1) cell_bounds = rect_coords( xmin, xmin + collen,