Make maximum rot vel a user input

documentation.md

@@ -248,10 +248,11 @@ The following fields are part of the floe settings (with default values):
  - min_floe_height: minimum floe height (0.1 m)
  - max_floe_height: maximum floe height (10.0 m)
  - min_aspect_ratio: minimum ratio between floe x-length and y-length by maximum coordiante values (0.05)
+ - maximum_ξ: the absolute maximum rotational velocity a floe can reach before it is capped at maximum_ξ (1e-5 rad/s)
  - nhistory: number of elements to save in floe's stress history (100)
  - subfloe_point_generator: generates floe's subfloe points (`MonteCarloPointsGenerator()`)
 
-If any of the minimum / maximum values are exceeded, a floe is removed in the course of the simulation.
+If any of the minimum values are exceeded, a floe is removed in the course of the simulation. If any of the maximum values are reached, the value is capped at the given value.
 
 There are two types of subfloe point generators.
 

src/physical_processes/process_settings.jl

@@ -26,6 +26,7 @@ Settings needed to create floes within the model.
  min_floe_height::FT = 0.1
  max_floe_height::FT = 10.0
  min_aspect_ratio::FT = 0.05
+ maximum_ξ::FT = 1e-5
  nhistory::Int = 100
  subfloe_point_generator::GT = MonteCarloPointsGenerator()
 
@@ -35,37 +36,36 @@ Settings needed to create floes within the model.
  min_floe_height,
  max_floe_height,
  min_aspect_ratio,
+ maximum_ξ,
  nhistory,
  subfloe_point_generator,
  ) where {FT <: AbstractFloat, GT <: AbstractSubFloePointsGenerator{FT}}
  if ρi < 0
- @warn "Ice density can't be negative. Resetting to default values \
- of 920."
+ @warn "Ice density can't be negative. Resetting to default values of 920."
  ρi = FT(920)
  end
  if min_floe_area < 0
- @warn "Floe area can't be negative. Resetting minimum floe area to \
- 0 m^2."
+ @warn "Floe area can't be negative. Resetting minimum floe area to 0 m^2."
  min_floe_area = FT(0)
  end
  if min_floe_height < 0
- @warn "Floe height can't be negative. Resetting minimum floe area \
- to 0,."
+ @warn "Floe height can't be negative. Resetting minimum floe area to 0."
  min_floe_height = FT(0)
  end
  if max_floe_height < 0
- @warn "Floe height can't be negative. Resetting maximum floe area \
- to default 10m."
+ @warn "Floe height can't be negative. Resetting to default of 10m."
  min_floe_height = FT(0)
  end
  if min_aspect_ratio < 0 || min_aspect_ratio > 1
- @warn "Aspect ratio must be between 0 and 1. Resetting to default \
- 0f 0.05."
+ @warn "Aspect ratio must be between 0 and 1. Resetting to default of 0.05."
+ min_aspect_ratio = FT(0.05)
+ end
+ if maximum_ξ < 0
+ @warn "Maximum rotational velocity must be greater than 0. Resetting to default of 1e-5."
  min_aspect_ratio = FT(0.05)
  end
  if nhistory < 1
- @warn "Need to save at least one timestep of stress history.
- Resetting to default of 100."
+ @warn "Need to save at least one timestep of stress history. Resetting to default of 100."
  nhistory = 100
  end
  new{FT, GT}(
@@ -74,6 +74,7 @@ Settings needed to create floes within the model.
  min_floe_height,
  max_floe_height,
  min_aspect_ratio,
+ maximum_ξ,
  nhistory,
  subfloe_point_generator,
  )
@@ -85,6 +86,7 @@ Settings needed to create floes within the model.
  min_floe_height,
  max_floe_height,
  min_aspect_ratio,
+ maximum_ξ,
  nhistory,
  subfloe_point_generator::GT,
  ) where {GT <: AbstractSubFloePointsGenerator} = 
@@ -94,6 +96,7 @@ Settings needed to create floes within the model.
  min_floe_height,
  max_floe_height,
  min_aspect_ratio,
+ maximum_ξ,
  nhistory,
  subfloe_point_generator,
  )

src/physical_processes/update_floe.jl

@@ -458,17 +458,17 @@ tendencies calculated at previous timesteps. Height, mass, stress, and strain
 also updated based on previous timestep thermodynamics and interactions with
 other floes. 
 Input:
- floe <Floe>
- Δt <Int> simulation timestep in second
- max_height <AbstractFloat> maximum floe height
+ floe  <Floe>
+ Δt  <Int> simulation timestep in second
+ floe_settings <FloeSettings> floe creation settings
 Output:
  None. Floe's fields are updated with values.
 """
 function timestep_floe_properties!(
  floes,
  tstep,
  Δt,
- max_height,
+ floe_settings,
 )
  Threads.@threads for i in eachindex(floes)
  cforce = floes.collision_force[i]
@@ -478,9 +478,9 @@ function timestep_floe_properties!(
  calc_stress!(LazyRow(floes, i))
  end
  # Ensure no extreem values due to model instability
- if floes.height[i] > max_height
- @warn "Reducing height to 10 m"
- floes.height[i] = max_height
+ if floes.height[i] > floe_settings.max_floe_height
+ @warn "Reducing height to $(floe_settings.max_floe_height) m"
+ floes.height[i] = floe_settings.max_floe_height
  end
 
  while maximum(abs.(cforce)) > floes.mass[i]/(5Δt)
@@ -536,9 +536,9 @@ function timestep_floe_properties!(
  end
  if frac != 1
  @warn "Adjusting u and v velocities to prevent too high"
+ dudt = frac*dudt
+ dvdt = frac*dvdt
  end
- dudt = frac*dudt
- dvdt = frac*dvdt
  floes.u[i] += 1.5Δt*dudt-0.5Δt*floes.p_dudt[i]
  floes.v[i] += 1.5Δt*dvdt-0.5Δt*floes.p_dvdt[i]
  floes.p_dudt[i] = dudt
@@ -547,9 +547,9 @@ function timestep_floe_properties!(
  dξdt = (floes.trqOA[i] + ctrq)/floes.moment[i]
  dξdt = frac*dξdt
  ξ = floes.ξ[i] + 1.5Δt*dξdt-0.5Δt*floes.p_dξdt[i]
- if abs(ξ) > 1e-5
+ if abs(ξ) > floe_settings.maximum_ξ
  @warn "Shrinking ξ" tstep = tstep
- ξ = sign(ξ) * 1e-5
+ ξ = sign(ξ) * floe_settings.maximum_ξ
  end
  floes.ξ[i] = ξ
  floes.p_dξdt[i] = dξdt

src/simulation_components/simulation.jl

@@ -165,7 +165,7 @@ function timestep_sim!(sim, tstep, start_tstep = 0)
  sim.model.floes,
  tstep,
  sim.Δt,
- sim.floe_settings.max_floe_height,
+ sim.floe_settings,
  )
  # Fracture floes
  if sim.fracture_settings.fractures_on && mod(tstep, sim.fracture_settings.Δt) == 0