diff --git a/documentation.md b/documentation.md index d1ee4dd..2a823dc 100644 --- a/documentation.md +++ b/documentation.md @@ -103,9 +103,9 @@ atmos = Atmos( ``` ### Domain: -The domain includes both the boundaries of the simulation and the topographic features. It defines the areas where the floes cannot go. Before you can define a domain, you need to define each boundary wall. As of now, only rectangular boundaries around the edge of the grid are allowed so we need a north, east, south, and west boundary wall. We have four types of boundary wall: open, collision, periodic, and compression. +The domain includes both the boundaries of the simulation and the topographic features. It defines the areas where the floes cannot go. Before you can define a domain, you need to define each boundary wall. As of now, only rectangular boundaries around the edge of the grid are allowed so we need a north, east, south, and west boundary wall. We have four types of boundary wall: open, collision, periodic, and moving. -With an open wall, if a floe overlaps with the boundary wall at all, the floe is removed from the simulation. With a collision wall, floes collide with the wall, and it is calculated similarly to a collision with another floe, except that since the wall cannot break, or be moved, it has an idealized force factor. With a periodic wall, if a floe overlaps with the wall, a copy of that floe is created passing back into the domain through the opposite wall. These two floes are equivalent and are linked. We call the copy a “ghost floe.” So, for example, if a floe is pushed partially out of the domain through the south wall by a current, a copy of the floe will be created, re-entering through the north wall. If one wall in the domain is periodic, its opposite wall must also be periodic, that is north-south and east-west. Compression boundaries are walls that can move with a constant velocity in either the x-direction (East and West walls) or in the y-direction (North and South walls). For example, if the west wall is a compression wall with a velocity of 0.1m/s it will move with a 0.1m/s u velocity and a 0m/s v velocity. These types of walls are for increasing pressure on the ice to investigate stress and strain on the floes. +With an open wall, if a floe overlaps with the boundary wall at all, the floe is removed from the simulation. With a collision wall, floes collide with the wall, and it is calculated similarly to a collision with another floe, except that since the wall cannot break, or be moved, it has an idealized force factor. With a periodic wall, if a floe overlaps with the wall, a copy of that floe is created passing back into the domain through the opposite wall. These two floes are equivalent and are linked. We call the copy a “ghost floe.” So, for example, if a floe is pushed partially out of the domain through the south wall by a current, a copy of the floe will be created, re-entering through the north wall. If one wall in the domain is periodic, its opposite wall must also be periodic, that is north-south and east-west. Moving boundaries are walls that can move with a constant velocity in either the x or y-direction, causing either a shear or compressive stress. For example, if the west wall is a moving wall with a x-velocity of 0.1m/s it will move with a 0.1m/s u velocity and a 0m/s v velocity towards the center of the domain. These types of walls are for increasing pressure on the ice to investigate stress and strain on the floes. Here is an example of creating a set of boundary walls using the grid from above: @@ -122,9 +122,9 @@ Once we have defined our walls, we can create a domain as follows: domain = Domain(nboundary, sboundary, eboundary, wboundary) ``` -Here is an example of how to create a compression wall as well that moves towards the center of the domain at 0.1m/s: +Here is an example of how to create a moving wall as well that moves towards the center of the domain with u velocity of 0m/s and a v velocity of -0.1m/s: ```julia -comp_boundary = CompressionBoundary(North, grid, -0.1) +move_boundary = MovingBoundary(North, grid, 0.0, -0.1) ``` diff --git a/examples/packed_compression.jl b/examples/moving_bounds.jl similarity index 84% rename from examples/packed_compression.jl rename to examples/moving_bounds.jl index c1406c1..3843bd5 100644 --- a/examples/packed_compression.jl +++ b/examples/moving_bounds.jl @@ -16,14 +16,14 @@ grid = RegRectilinearGrid( Δgrid, Δgrid, ) -ocean = Ocean(grid, 0.1, 0.0, 0.0) +ocean = Ocean(grid, 0.0, 0.0, 0.0) atmos = Atmos(grid, 0.0, 0.0, -1.0) # Domain creation -nboundary = CompressionBoundary(North, grid, -0.1) -sboundary = CompressionBoundary(South, grid, 0.1) -eboundary = CollisionBoundary(East, grid) -wboundary = CollisionBoundary(West, grid) +nboundary = MovingBoundary(North, grid, 0.2, 0.0) +sboundary = MovingBoundary(South, grid, 0.2, 0.0) +eboundary = PeriodicBoundary(East, grid) +wboundary = PeriodicBoundary(West, grid) domain = Domain(nboundary, sboundary, eboundary, wboundary) @@ -37,13 +37,15 @@ floe_arr = initialize_floe_field( Δh, rng = Xoshiro(1), ) - +nfloes = length(floe_arr) +floe_arr.u .= 0 +floe_arr.v .= -0.01 # 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) +consts = Constants(E = modulus, Cd_io = 0.0, f = 0.0, turnθ = 0.0) # Run simulation run_time!(simulation) = @time run!(simulation) @@ -80,12 +82,11 @@ simulation = Simulation( model = model, consts = consts, Δt = Δt, - nΔt = 10000, + nΔt = 5000, verbose = true, writers = writers, rng = Xoshiro(1), coupling_settings = coupling_settings, - weld_settings = weld_settings, ) run_time!(simulation) diff --git a/src/Subzero.jl b/src/Subzero.jl index 06c797d..cf377e4 100644 --- a/src/Subzero.jl +++ b/src/Subzero.jl @@ -13,7 +13,7 @@ export CollisionBC, CompressionBC, AbstractBoundary, - CompressionBoundary, + MovingBoundary, CollisionBoundary, PeriodicBoundary, OpenBoundary, diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl index 3c26877..cf2d1d7 100644 --- a/src/physical_processes/collisions.jl +++ b/src/physical_processes/collisions.jl @@ -204,11 +204,17 @@ boundary. """ get_velocity( element::AbstractDomainElement{FT}, - x, - y, + _, + _, ) where {FT} = FT(0), FT(0) +get_velocity( + element::MovingBoundary, + _, + _, +) = element.u, element.v + """ calc_friction_forces( ifloe, @@ -603,7 +609,7 @@ function floe_domain_element_interaction!( floe, element::Union{ CollisionBoundary, - CompressionBoundary, + MovingBoundary, TopographyElement, }, elem_idx, @@ -676,17 +682,17 @@ end Move North/South compression boundaries by given velocity. Update coords and val fields to reflect new position. Inputs: - boundary + 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::CompressionBoundary{D, FT}, + boundary::MovingBoundary{D, FT}, Δt, ) where {D <: Union{North, South}, FT <: AbstractFloat} - Δd = boundary.velocity * Δt + Δd = boundary.v * Δt boundary.val += Δd translate!(boundary.coords, FT(0), Δd) end @@ -696,17 +702,17 @@ end Move East/West compression boundaries by given velocity. Update coords and val fields to reflect new position. Inputs: - boundary + boundary domain compression boundary Δt number of seconds in a timestep Outputs: None. Move boundary East/West depending on velocity. """ function update_boundary!( - boundary::CompressionBoundary{D, FT}, + boundary::MovingBoundary{D, FT}, Δt, ) where {D <: Union{East, West}, FT <: AbstractFloat} - Δd = boundary.velocity * Δt + Δd = boundary.u * Δt boundary.val += Δd translate!(boundary.coords, Δd, FT(0)) end diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index c068419..b892df5 100644 --- a/src/physical_processes/coupling.jl +++ b/src/physical_processes/coupling.jl @@ -608,7 +608,7 @@ function in_bounds( end """ - calc_mc_values!( + calc_subfloe_values!( floe::Union{Floe{FT}, LazyRow{Floe{FT}}}, grid, domain, @@ -616,7 +616,7 @@ end mc_grid_idx, ) -Calculates monte carlo point's cartesian coordiantes, polar coordiantes, +Calculates subfloe point's cartesian coordiantes, polar coordiantes, velocity and index within the grid. Inputs: floe floe @@ -635,7 +635,7 @@ Outputs: mc_cart and mc_grid_idx filled with data for given floe's monte carlo points up to row j. """ -function calc_mc_values!( +function calc_subfloe_values!( floe::Union{Floe{FT}, LazyRow{Floe{FT}}}, grid, domain, @@ -653,7 +653,7 @@ function calc_mc_values!( x = px + floe.centroid[1] # at centroid y = py + floe.centroid[2] # at centroid # If point is in bounds, continue to find rest of values - if in_bounds(x, y, grid, domain.east, domain.north) + if in_bounds(x, y, grid, domain.north, domain.east) j += 1 # if added to outputs, move to next index in output array cart_vals[j, 1] = x cart_vals[j, 2] = y @@ -1509,7 +1509,7 @@ function calc_one_way_coupling!( grid_idx = Matrix{Int}(undef, max_points, 2) for i in eachindex(floes) # Monte carlo point cartesian coordinates and grid cell indices - npoints = calc_mc_values!( + npoints = calc_subfloe_values!( LazyRow(floes, i), grid, domain, diff --git a/src/simulation_components/domain_and_grid.jl b/src/simulation_components/domain_and_grid.jl index 4206749..dc5c28c 100644 --- a/src/simulation_components/domain_and_grid.jl +++ b/src/simulation_components/domain_and_grid.jl @@ -338,66 +338,74 @@ function CollisionBoundary{D, FT}( end """ - CompressionBoundary <: AbstractBoundary + 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 velocity is in [m/s]. +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 CompressionBoundary{D, FT}<:AbstractBoundary{D, FT} +mutable struct MovingBoundary{D, FT}<:AbstractBoundary{D, FT} coords::PolyVec{FT} val::FT - velocity::FT + u::FT + v::FT end """ - CompressionBoundary(::Type{FT}, ::Type{D}, args...) + MovingBoundary(::Type{FT}, ::Type{D}, args...) -A float type FT can be provided as the first argument of any CompressionBoundary +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 CompressionBoundary of type FT and D will be created by passing all other +A MovingBoundary of type FT and D will be created by passing all other arguments to the correct constructor. """ -CompressionBoundary(::Type{FT}, ::Type{D}, args...) where { +MovingBoundary(::Type{FT}, ::Type{D}, args...) where { FT <: AbstractFloat, D <: AbstractDirection, -} = CompressionBoundary{D, FT}(args...) +} = MovingBoundary{D, FT}(args...) """ - CompressionBoundary(::Type{D}, args...) + MovingBoundary(::Type{D}, args...) -If a float type isn't specified, CompressionBoundary will be of type Float64 and +If a float type isn't specified, MovingBoundary will be of type Float64 and the correct constructor will be called with all other arguments. """ -CompressionBoundary(::Type{D}, args...) where {D <: AbstractDirection} = - CompressionBoundary{D, Float64}(args...) +MovingBoundary(::Type{D}, args...) where {D <: AbstractDirection} = + MovingBoundary{D, Float64}(args...) """ - CompressionBoundary{D, FT}(grid, velocity) + 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 - velocity velocity of boundary + u u velocity of boundary + v v velocity of boundary Outputs: - CompressionBoundary on edge of grid given by direction. -Note: - If a boundary should move from north to south or east to west, the velocity - should be negative. Else the velocity should be positive. + MovingBoundary on edge of grid given by direction. """ -function CompressionBoundary{D, FT}( +function MovingBoundary{D, FT}( grid, - velocity, + u, + v, ) where {D <: AbstractDirection, FT <: AbstractFloat} val, coords = boundary_coords(grid, D) - CompressionBoundary{D, FT}( + MovingBoundary{D, FT}( coords, val, - velocity, + u, + v, ) end @@ -409,7 +417,7 @@ Union of all non-peridic boundary types to use as shorthand for dispatch. const NonPeriodicBoundary = Union{ OpenBoundary, CollisionBoundary, - CompressionBoundary, + MovingBoundary, } """ diff --git a/test/test_model.jl b/test/test_model.jl index 5a33184..f7cde87 100644 --- a/test/test_model.jl +++ b/test/test_model.jl @@ -176,7 +176,7 @@ b2 = Subzero.OpenBoundary(East, g) b3 = Subzero.CollisionBoundary(West, g) b4 = Subzero.PeriodicBoundary(South, g) - b5 = Subzero.CompressionBoundary(South, g, 1.0) + 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]]] @@ -189,8 +189,8 @@ @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.velocity == 1.0 - @test b5.val == 0.0 + @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 diff --git a/test/test_physical_processes/test_collisions.jl b/test/test_physical_processes/test_collisions.jl index 7e5b18f..c18dad7 100644 --- a/test/test_physical_processes/test_collisions.jl +++ b/test/test_physical_processes/test_collisions.jl @@ -295,13 +295,13 @@ @test all(corner_floe.interactions[:, xforce] .<= 0) @test all(corner_floe.interactions[:, yforce] .<= 0) # Test compression boundaries movement - nc_boundary = CompressionBoundary(North, grid, -0.1) + nc_boundary = MovingBoundary(North, grid, 0.0, -0.1) nc_coords = deepcopy(nc_boundary.coords) - sc_boundary = CompressionBoundary(South, grid, 0.1) + sc_boundary = MovingBoundary(South, grid, 0.0, 0.1) sc_coords = deepcopy(sc_boundary.coords) - ec_boundary = CompressionBoundary(East, grid, 0.1) + ec_boundary = MovingBoundary(East, grid, 0.1, 0.0) ec_coords = deepcopy(ec_boundary.coords) - wc_boundary = CompressionBoundary(West, grid, 0.1) + 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)