Moving Boundaries (#86)

* Update compression bounds to moving bounds

* Cleanup comments

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)
 ```
 
 

examples/packed_compression.jl → 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)
 

src/Subzero.jl

@@ -13,7 +13,7 @@ export
  CollisionBC,
  CompressionBC,
  AbstractBoundary,
- CompressionBoundary,
+ MovingBoundary,
  CollisionBoundary,
  PeriodicBoundary,
  OpenBoundary,

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 <CompressionBoundary{Union{North, South}, AbstractFloat}> 
+ boundary <MovingBoundary{Union{North, South}, AbstractFloat}> 
  domain compression boundary
  Δt <Int> 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 <CompressionBoundary{Union{East, West}, AbstractFloat}> 
+ boundary <MovingBoundary{Union{East, West}, AbstractFloat}> 
  domain compression boundary
  Δt <Int> 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

src/physical_processes/coupling.jl

@@ -608,15 +608,15 @@ function in_bounds(
 end
 
 """
- calc_mc_values!(
+ calc_subfloe_values!(
  floe::Union{Floe{FT}, LazyRow{Floe{FT}}},
  grid,
  domain,
  mc_cart,
  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 <Union{Floe{AbstractFloat}, LazyRow{Floe{AbstractFloat}}}> 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,

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 <AbstractGrid> model grid
- velocity <AbstractFloat> velocity of boundary
+ u <AbstractFloat> u velocity of boundary
+ v <AbstractFloat> 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,
 }
 
 """

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

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)