Clarify stress calcultor comments and functions

src/physical_processes/fractures.jl

@@ -300,8 +300,6 @@ end
 # This can be changed to add some sort of scaling to the DamageStressCalculator, potentially
 # based on area like the function above.
 function scale_stress!(::DamageStressCalculator, σvals, min_floe_area, floe_area)
- @warn "Current implementation of DamageStressCalculator is unfinished. Currently, stress
- is not being scaled."
  return
 end 
 
@@ -401,7 +399,6 @@ function split_floe(
  floe_settings,
  Δt,
 ) where {FT}
- τ = floe_settings.stress_calculator.τ
  new_floes = StructArray{Floe{FT}}(undef, 0)
  # Generate voronoi tesselation in floe's bounding box
  scale_fac = fill(2floe.rmax, 2)

src/physical_processes/types.jl

@@ -12,18 +12,31 @@ abstract type AbstractStressCalculator{FT<:AbstractFloat} end
 """
  DecayAreaScaledCalculator{FT<:AbstractFloat}<:AbstractStressCalculator
 
-Type of AbstractStressCalculator that implements stress calculations the same way that 
-Brandon Montemuro and Georgy Manucharyan do in the MatLab versio of the code. The area-scaling
-part of the stress calculations keep a record of previous stress (and implicity damage). The
-decay aspect increases importance placed on new damage.
+Type of AbstractStressCalculator that implements stress calculations by accumulating each
+timestep of stress using a decay equation. The decay aspect increases importance placed on
+new damage. The decay equation is as follows:
+
+stress_accum = stress_accum + Δt/τ .* (stress_instant - stress_accum)
+
+To avoid recomputation, Δt/τ is saved in the floe's damage field. 
+
+The area-scaling part of the stress calculations comes into play when deciding if a floe
+will be fractured. The accumulated stress can be scaled by a ratio of
+(floe.area/min_floe_area)^α. By changing the α value, either larger or smaller floes
+fracture more eassily. The scaled value will not be saved as this is equivalent to scaling
+the simulation fracture criteria (morh's cone or hibler's ellipse, etc.), but it is less
+computationally intensive. 
+
 Fields:
- τ <AbstractFloat> Difference between current and previous stress scaled by Δt/τ
+ λ <AbstractFloat> Decay parameter used when calculating accumulated stress in
+ update_stress! function. Commonly set to Δt/τ. Should be between 0 - 1 such
+ that stress_accum = stress_accum(1-λ) + stress_instant(λ)
  α <AbstractFloat> Adjusts ellipse in stress space by raising area ratio to the α
 Note:
- τ is used in calc_stress!(), whereas α is used in determine_fractures().
+ τ is used in calc_stress!, whereas α is used in determine_fractures.
 """
 @kwdef struct DecayAreaScaledCalculator{FT<:AbstractFloat} <: AbstractStressCalculator{FT}
- τ::FT = 25.0
+ λ::FT = 0.2
  α::FT = 0.0
 end
 
@@ -52,24 +65,35 @@ DecayAreaScaledCalculator(args...; kwargs...) = DecayAreaScaledCalculator{Float6
 """
  DamageStressCalculator
 
-Type of AbstractStressCalculator that calculates stress with damage*currstress, as suggested
-by Mukund Gupta. This method keeps track of damage directly and rather than modifying the 
-"accumulated stress" to store damage, as done in DecayAreaScaledCalculator. This calculator
-keeps track of an explicit parameter for damage.
+Type of AbstractStressCalculator that calculates stress with damage * stress_instant, as
+suggested by Mukund Gupta. This method could keep track of damage directly, perhapes as a
+value between 0-1, and rather than calculating an "accumulated stress" to store damage, as
+done in DecayAreaScaledCalculator.
+
+In this calculator, the floe's damage field keeps track of an explicit parameter for damage.
 Fields:
  τ <AbstractFloat> Difference between current and previous stress scaled by Δt/τ.
  This field should represent damage.
 Note:
  This method is not fully implemented. The infrastructure for the damage calculator is
- provided, but functions that depend on this calculator,  
+ provided, but functions that depend on this calculator need to be implemented. 
  There is no α parameter in this calculator because currently it does not adjust the
  boundary in stress space by multiplying the eigenvalues of stress_accum by something.
  This could be implemented if the user desires.
+
+ The functions that need to be implemented are: scale_stress!, update_damage!, and
+ update_stress!. You may also need to create an inialize damage parameter, which can set
+ an initial damage based on the calculator type if you don't want it to start at 0.
+
+ You can add needed fields to the below struct as in the DecayAreaScaledCalculator
 """
 @kwdef struct DamageStressCalculator{FT<:AbstractFloat} <: AbstractStressCalculator{FT}
- τ::FT = 20.0
-end
+ function DamageStressCalculator{FT}() where FT
+ throw(ErrorException("DamageStressCalculator isn't implemented yet!! Take a look at the issue/comments if you're interested in implementing."))
+ return new{FT}()
+ end
 
+end
 
 """
  DamageStressCalculator(::Type{FT}; kwargs...)

src/physical_processes/update_floe.jl

@@ -394,51 +394,47 @@ function calc_stress!(floe::FloeType{FT}, floe_settings, Δt) where {FT}
  xi, yi = floe.centroid
  inters = floe.interactions
  # Calculates timestep stress
- stress = fill(FT(0), 2, 2)
- for i in 1:floe.num_inters
- stress[1, 1] += (inters[i, xpoint] - xi) * inters[i, xforce]
- stress[1, 2] += (inters[i, ypoint] - yi) * inters[i, xforce] +
- (inters[i, xpoint] - xi) * inters[i, yforce]
- stress[2, 2] += (inters[i, ypoint] - yi) * inters[i, yforce]
+ stress = zeros(FT, 2, 2)
+ if floe.num_inters > 0
+ for i in 1:floe.num_inters
+ stress[1, 1] += (inters[i, xpoint] - xi) * inters[i, xforce]
+ stress[1, 2] += (inters[i, ypoint] - yi) * inters[i, xforce] +
+ (inters[i, xpoint] - xi) * inters[i, yforce]
+ stress[2, 2] += (inters[i, ypoint] - yi) * inters[i, yforce]
+ end
+ stress[1, 2] *= FT(0.5)
+ stress[2, 1] = stress[1, 2]
+ stress .*= 1/(floe.area * floe.height)
  end
- stress[1, 2] *= FT(0.5)
- stress[2, 1] = stress[1, 2]
- stress .*= 1/(floe.area * floe.height)
- update_damage!(floe_settings.stress_calculator, stress, floe, Δt)
+ update_damage!(floe_settings.stress_calculator, stress, floe)
  update_stress!(floe_settings.stress_calculator, stress, floe)
- floe.stress_instant = stress
  return
 end
 
-# This function updates the damage parameter based on the type of calculator we are working
-# with. It is empty because the DamageStressCalculator is not fully implemented
-function update_damage!(stress_calculator, curr_stress, floe, Δt)
+# The damage parameter isn't used with the DecayAreaScaledCalculator.
+function update_damage!(stress_calculator::DecayAreaScaledCalculator, curr_stress, floe)
  return
 end
 
-# This is where one would implement a different method of updating the damage parameter.
-# currently this effectively carries out the same math as update_stress! for DecayAreaScaledCalculator
-function update_damage!(stress_calculator::DamageStressCalculator, curr_stress, floe, Δt)
- @warn "Current implementation of DamageStressCalculator is unfinished. Calculating 
- stress as if using DecayAreaScaledCalculator"
- τ = stress_calculator.τ
- if iszero(curr_stress)
- floe.damage = 0
- else
- floe.damage = (floe.stress_accum + (Δt/τ)*(curr_stress - floe.stress_accum))./curr_stress
- end
+#= This is where one would implement a different method of updating the damage parameter.
+It is empty because the DamageStressCalculator is not fully implemented but it should
+update the floe.damage parameter. =#
+function update_damage!(stress_calculator::DamageStressCalculator, curr_stress, floe)
+ return
 end
 
-# Updating stress according to DamageStressCalculator. When DamageStressCalculator is 
-# fully implemented, this function does not have to change.
-function update_stress!(stress_calculator::DamageStressCalculator, curr_stress, floe)
- floe.stress_accum = floe.damage .* curr_stress
+# Updating stress using a decay equation where the curr_stress is λ
+function update_stress!(stress_calculator::DecayAreaScaledCalculator, curr_stress, floe)
+ λ = stress_calculator.λ
+ floe.stress_accum = (1 - λ) * floe.stress_accum + λ * curr_stress
+ floe.stress_instant = curr_stress
+ return
 end
 
-# Updating stress in the same way that Brandon Montemuro and Georgy Manucharyan do 
-# in their MatLab code.
-function update_stress!(stress_calculator, curr_stress, floe)
- floe.stress_accum = floe.stress_accum + floe.damage .* (curr_stress - floe.stress_accum)
+#= Updating stress according to DamageStressCalculator. This function should use the
+floe.damage parameter. =#
+function update_stress!(stress_calculator::DamageStressCalculator, curr_stress, floe)
+ return
 end
 
 """

src/simulation_components/floe.jl

@@ -66,7 +66,7 @@ Singular sea ice floe with fields describing current state.
  stress_accum::Matrix{FT} = zeros(2, 2)
  stress_instant::Matrix{FT} = zeros(2, 2)
  strain::Matrix{FT} = zeros(2, 2)
- damage::Matrix{FT} = zeros(2, 2)
+ damage::FT = 0.0 # damage to the floe (to be used in new stress calculator)
  # Previous values for timestepping -------------------------------------
  p_dxdt::FT = 0.0 # previous timestep x-velocity
  p_dydt::FT = 0.0 # previous timestep y-velocity
@@ -293,7 +293,6 @@ function poly_to_floes!(
  poly::Polys,
  hmean,
  Δh;
- damage = calc_initial_damage(floe_settings.stress_calculator, Δt),
  floe_settings = floe_settings,
  rng = rng,
  kwargs...
@@ -314,21 +313,6 @@ function poly_to_floes!(
  return 0
 end
 
-# When using DecayAreaScaledCalculator, we initialize stress to dt/τ, which is the timestep
-# over a parameter that we set when creating the stress calculator. This setup is 
-# necessary to calculate stress in the same way that Georgy Manucharyan and Brandon Montemuro
-# do in their code.
-calc_initial_damage(stress_calculator::DecayAreaScaledCalculator, Δt) = fill(Δt / stress_calculator.τ, 2, 2)
-
-# This is a place holder untio DamageStressCalculator is fully implemented. This should be
-# replaced with a physical interpretation of what the damage parameter should start as.
-# Damage does not necessary need to be a matrix, it has just been implemented this way to
-# be able to mimic Manucharyan's and Montemuro's code.
-function calc_initial_damage(stress_calculator::DamageStressCalculator, Δt)
- @warn "DamageStressCalculator not fully implemented. Initializing damage to zero." 
- return zeros(2, 2)
-end
-
 """
  initialize_floe_field(args...)
 
@@ -571,7 +555,6 @@ function initialize_floe_field(
  floe_settings = FloeSettings(FT, min_floe_area = 0),
  rng = Xoshiro(),
 ) where {FT <: AbstractFloat}
- τ = floe_settings.stress_calculator.τ
  floe_arr = StructArray{Floe{FT}}(undef, 0)
  nfloes_added = 0
  # Availible space in domain