Add floe settings

README.md

@@ -5,7 +5,7 @@
 
 <!-- description -->
 <p align="center">
- <strong>:ice_cube: Sea-ice model to explicitly simulate individual floe life cycles using complex discrete elements with time-evolving shapes. Easy to couple with Oceananigans for explorations of coupled ice-ocean dynamics. Ported from MATLAB to Julia. :ice_cube:</strong>
+ <strong> Sea-ice model to explicitly simulate individual floe life cycles using complex discrete elements with time-evolving shapes. Easy to couple with Oceananigans for explorations of coupled ice-ocean dynamics. Ported from MATLAB to Julia.</strong>
 </p>
 
 <!-- Information badges -->
@@ -63,7 +63,7 @@ which will add Subzero as a package.
 
 ## Running your first model:
 
-For detailed instructions on how to create different types of models and simulations, please see the [documentation.md file](https://github.com/Caltech-OCTO/Subzero.jl/blob/documentation/documentation.md). However, we will give a basic example here.
+For detailed instructions on how to create different types of models and simulations, please see the [documentation.md file](https://github.com/Caltech-OCTO/Subzero.jl/blob/main/documentation.md). However, we will give a basic example here.
 
 Let’s run a basic simulation with initially stationary floes pushed into a collision boundary by a uniform, zonally flowing ocean. In this simulation, collisions between floes are on by default and we will enable floe fracturing. 
 
@@ -85,59 +85,74 @@ domain = Domain(
  CollisionBoundary(West, grid),
 )
 # Create floes
+floe_settings = FloeSettings(nhistory = 100, min_floe_area = 1e5)
 floe_arr = initialize_floe_field(
- 0, # number of floes
+ Float64,
+ 100, # number of floes
  [0.7], # floe concentration
  domain,
  0.5, # average floe height
- 0.05, # floe height variability
+ 0.05; # floe height variability
+ floe_settings = floe_settings,
 ) 
 # Create model
 model = Model(grid, ocean, atmos, domain, floe_arr) 
-# Create simulation
+
+# Create settings
 modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area))) 
 consts = Constants(E = modulus) 
+
 fracture_settings = FractureSettings( 
  fractures_on = true,
  criteria = HiblerYieldCurve(floe_arr),
  Δt = 75,
  npieces = 3,
- nhistory = 1000,
  deform_on = true, 
 )
-
+# Create output writers
+dir = "output/sim"
 initwriter = InitialStateOutputWriter(
  dir = dir,
  filename = "initial_state.jld2",
  overwrite = true,
  )
 floewriter = FloeOutputWriter(
  100,
- dir = "output/sim",
+ dir = dir,
  filename = "floes.jld2",
  overwrite = true,
 )
 writers = OutputWriters(initwriter, floewriter)
+# Create simulation
+Δt = 10
 simulation = Simulation( 
  model = model, 
  consts = consts, 
- Δt = Δt, 
- nΔt = 5000, 
- verbose = true, 
- fracture_settings,
+ Δt = Δt,  # 10 second timesteps 
+ nΔt = 10000, # Run for 10000 timesteps
+ verbose = true; 
+ fracture_settings = fracture_settings,
  writers = writers,
 )
 # Run simulation
 run!(simulation)
+
+# Plot simulation
+plot_sim(
+ "output/sim/floes.jld2",
+ "output/sim/initial_state.jld2",
+ Δt,
+ "output/sim/example.mp4",
+)
 ``` 
 
-Check out our [documentation](https://github.com/Caltech-OCTO/Subzero.jl/blob/documentation/documentation.md) for more examples and explanations for the code above.
+Check out our [documentation](https://github.com/Caltech-OCTO/Subzero.jl/blob/main/documentation.md) for more examples and explanations for the code above.
 
 ## Citing:
 We still need to figure this out. Please reach out so we can discuss.
 
 ## Contributing:
-If you’re interested in helping develop Subzero, have found a bug, or have a new feature that you want implemented, please open an issue on the repository and we can talk about this. We will be working on a contributers' guide for this in the future.
+If you’re interested in helping develop Subzero, have found a bug, or have a new feature that you want implemented, please open an issue on the repository and we can talk about this. We will be working on a contributers' guide in the future.
 
 ## Movies:
 

examples/shear_flow.jl

@@ -36,9 +36,11 @@ wboundary = PeriodicBoundary(West, grid)
 domain = Domain(nboundary, sboundary, eboundary, wboundary)
 
 coupling_settings = CouplingSettings(
- subfloe_point_generator = SubGridPointsGenerator(grid, 2),
  two_way_coupling_on = true,
 )
+floe_settings = FloeSettings(
+ subfloe_point_generator = SubGridPointsGenerator(grid, 2),
+)
 # Floe creation
 floe_arr = initialize_floe_field(
  FT,
@@ -48,7 +50,7 @@ floe_arr = initialize_floe_field(
  hmean,
  Δh;
  rng = Xoshiro(1),
- coupling_settings = coupling_settings
+ floe_settings = floe_settings
 )
 
 # Model creation
@@ -70,7 +72,7 @@ simulation = Simulation(
  model = model,
  consts = consts,
  Δt = Δt,
- nΔt = 1000,
+ nΔt = 10000,
  verbose = true,
  writers = writers,
  rng = Xoshiro(1),

examples/simple_strait.jl

@@ -37,9 +37,11 @@ topo_arr = initialize_topography_field(
 domain = Domain(nboundary, sboundary, eboundary, wboundary, topo_arr)
 
 coupling_settings = CouplingSettings(
- subfloe_point_generator = SubGridPointsGenerator(grid, 2),
  two_way_coupling_on = true,
 )
+floe_settings = FloeSettings(
+ subfloe_point_generator = SubGridPointsGenerator(grid, 2),
+)
 
 # Floe creation
 floe_arr = initialize_floe_field(
@@ -50,15 +52,14 @@ floe_arr = initialize_floe_field(
  hmean,
  Δh,
  rng = Xoshiro(3),
- coupling_settings = coupling_settings,
+ floe_settings = floe_settings,
 )
 
 fracture_settings = FractureSettings(
  fractures_on = true,
  criteria = HiblerYieldCurve(floe_arr),
  Δt = 75,
  npieces = 3,
- nhistory = 1000,
  deform_on = false,
 )
 
@@ -87,6 +88,7 @@ simulation = Simulation(
  Δt = Δt,
  nΔt = 10000,
  verbose = true,
+ floe_settings = floe_settings,
  coupling_settings = coupling_settings,
  fracture_settings = fracture_settings,
  ridgeraft_settings = ridgeraft_settings,

examples/test_run.jl

@@ -153,8 +153,10 @@ coords = [[
  [0.0, 0.0],
 ]]
 
-coupling_settings = CouplingSettings(
+floe_settings = FloeSettings(
  subfloe_point_generator = SubGridPointsGenerator(grid, 1),
+)
+coupling_settings = CouplingSettings(
  two_way_coupling_on = true,
 )
 
@@ -165,7 +167,7 @@ floe_arr = initialize_floe_field(
  domain,
  0.5,
  0.0;
- coupling_settings = coupling_settings,
+ floe_settings = floe_settings,
  rng = Xoshiro(1),
 )
 floe_arr.α .= 1e-5
@@ -197,6 +199,7 @@ simulation = Simulation(
  verbose = true,
  consts = Constants(Cd_ia = 0, Cd_ao = 0),
  coupling_settings = coupling_settings,
+ floe_settings = floe_settings,
 
 )
 

src/Subzero.jl

@@ -58,6 +58,7 @@ export
  SimplificationSettings,
  RidgeRaftSettings,
  WeldSettings,
+ FloeSettings,
  PolyVec,
  OutputWriters,
  check_energy_momentum_conservation_julia,

src/physical_processes/fractures.jl

@@ -334,8 +334,7 @@ function deform_floe!(
  floe,
  deformer_coords,
  deforming_forces,
- coupling_settings,
- consts,
+ floe_settings,
  Δt,
  rng,
 )
@@ -366,8 +365,7 @@ function deform_floe!(
  floe,
  new_floe_poly,
  floe.mass,
- coupling_settings,
- consts,
+ floe_settings,
  rng,
  )
  conserve_momentum_change_floe_shape!(
@@ -388,8 +386,7 @@ end
  floe,
  rng,
  fracture_settings,
- coupling_settings,
- consts,
+ floe_settings,
  Δt,
  )
 Splits a given floe into pieces using voronoi tesselation.
@@ -398,8 +395,7 @@ Inputs:
  floe <Floe> floe in simulation
  rng <RNG> random number generator used for voronoi tesselation
  fracture_settings <FractureSettings> simulation's fracture settings
- coupling_settings <CouplingSettings> simulation's coupling settings
- consts <Constants> simulation's constants
+ floe_settings <FloeSettings> simulation's settings for making floes
  Δt <Int> length of simulation timesteps in seconds
 Outputs:
  new_floes <StructArray{Floes}> list of pieces floe is split into, each of
@@ -409,8 +405,7 @@ function split_floe(
  floe::Union{Floe{FT}, LazyRow{Floe{FT}}},
  rng,
  fracture_settings,
- coupling_settings,
- consts,
+ floe_settings,
  Δt,
 ) where {FT}
  new_floes = StructArray{Floe{FT}}(undef, 0)
@@ -444,15 +439,13 @@ function split_floe(
  for i in eachindex(pieces_polys)
  if pieces_areas[i] > 0
  mass = floe.mass * (pieces_areas[i]/total_area)
- height = mass / (consts.ρi * pieces_areas[i])
+ height = mass / (floe_settings.ρi * pieces_areas[i])
  pieces_floes = poly_to_floes(
  FT,
  pieces_polys[i],
  height,
  0; # Δh - range of random height difference between floes
- ρi = consts.ρi,
- coupling_settings = coupling_settings,
- fracture_settings = fracture_settings,
+ floe_settings = floe_settings,
  rng = rng,
  u = floe.u,
  v = floe.v,
@@ -482,9 +475,7 @@ end
  max_floe_id,
  rng,
  fracture_settings,
- coupling_settings,
- simp_settings,
- consts::Constants{FT},
+ floe_settings,
  Δt,
  )
 Fractures floes that meet the criteria defined in the fracture settings.
@@ -493,9 +484,7 @@ Inputs:
  max_floe_id <Int> maximum ID of any floe created so far in simulation
  rng <RNG> random number generator
  fracture_settings <FractureSettings> sim's fracture settings
- coupling_settings <CouplingSettings> sim's coupling settings
- simp_settings <SimplificationSettings> sim's simplification settings
- consts <Constants> sim's constants
+ floe_settings <FloeSettings> sim's settings to make floes
  Δtout <Int> length of simulation timestep in seconds
 Outputs:
  max_floe_id <Int> new highest floe ID after adding new floes to floe array.
@@ -506,16 +495,14 @@ function fracture_floes!(
  max_floe_id,
  rng,
  fracture_settings,
- coupling_settings,
- simp_settings,
- consts,
+ floe_settings,
  Δt,
 ) where {FT <: AbstractFloat}
  # Determine which floes will fracture
  frac_idx = determine_fractures(
  floes,
  fracture_settings.criteria,
- simp_settings.min_floe_area,
+ floe_settings.min_floe_area,
  )
  # Initialize list for new floes created from fracturing existing floes
  nfloes2frac = length(frac_idx)
@@ -544,8 +531,7 @@ function fracture_floes!(
  LazyRow(floes, frac_idx[i]), 
  floes.coords[deforming_floe_idx],
  deforming_inter[xforce:yforce],
- coupling_settings,
- consts,
+ floe_settings,
  Δt,
  rng,
  )
@@ -557,8 +543,7 @@ function fracture_floes!(
  LazyRow(floes, frac_idx[i]),
  rng,
  fracture_settings,
- coupling_settings,
- consts,
+ floe_settings,
  Δt,
  )
  append!(fracture_list[i], new_floes)