From 3402de71508fedbd959f8287e7b375f21b67171c Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Mon, 17 Jul 2023 17:35:29 -0700 Subject: [PATCH] Debug type inference (#76) * Start debugging type instability * Reduce type instability * Attempt to update ghost floes * Fix ghost bug * Work on removing deepcopy * Remove output type instabilities * Remove direct floe modifications * Fix floe util * Update intersect coords * Remove stuct array i index call * Update comments * Finish tests --- Project.toml | 4 +- examples/shear_flow.jl | 6 +- examples/simple_strait.jl | 5 +- examples/test_run.jl | 106 ++-- src/floe_utils.jl | 212 ++++--- src/output.jl | 99 ++-- src/physical_processes/collisions.jl | 520 +++++++++++------- src/physical_processes/coupling.jl | 13 +- src/physical_processes/update_floe.jl | 9 +- src/simulation_components/domain_and_grid.jl | 19 +- src/simulation_components/model.jl | 21 +- src/simulation_components/simulation.jl | 20 +- test/test_floe_utils.jl | 15 +- .../test_collisions.jl | 1 + 14 files changed, 632 insertions(+), 418 deletions(-) diff --git a/Project.toml b/Project.toml index 2e10a9a..17c5c26 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" @@ -14,14 +15,15 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PolygonInbounds = "e4521ec6-8c1d-418e-9da2-b3bc4ae105d6" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/examples/shear_flow.jl b/examples/shear_flow.jl index 6ba7fd9..098c14b 100644 --- a/examples/shear_flow.jl +++ b/examples/shear_flow.jl @@ -1,4 +1,4 @@ -using JLD2, Random, Statistics, StructArrays, Subzero +using JLD2, Random, Statistics, StructArrays, Subzero, ProfileView import LibGEOS as LG # User Inputs @@ -71,13 +71,15 @@ simulation = Simulation( model = model, consts = consts, Δt = Δt, - nΔt = 2000, + nΔt = 1000, verbose = true, writers = writers, rng = Xoshiro(1), coupling_settings = coupling_settings, ) run_time!(simulation) + +#ProfileView.@profview run!(simulation) Subzero.create_sim_gif("output/shear_flow/floes.jld2", "output/shear_flow/initial_state.jld2", diff --git a/examples/simple_strait.jl b/examples/simple_strait.jl index 8cdb863..91e9481 100644 --- a/examples/simple_strait.jl +++ b/examples/simple_strait.jl @@ -1,11 +1,11 @@ -using JLD2, Random, SplitApplyCombine, Statistics, StructArrays, Subzero +using JLD2, Random, SplitApplyCombine, Statistics, StructArrays, Subzero, ProfileView import LibGEOS as LG # User Inputs const FT = Float64 const Lx = 1e5 const Ly = 1e5 -const Δgrid = 1e4 +const Δgrid = 2e3 const hmean = 0.25 const Δh = 0.0 const Δt = 20 @@ -90,6 +90,7 @@ simulation = Simulation( # Run simulation run_time!(simulation) +#ProfileView.@profview run!(simulation) Subzero.create_sim_gif("output/simple_strait/floes.jld2", "output/simple_strait/initial_state.jld2", diff --git a/examples/test_run.jl b/examples/test_run.jl index 3f4fc62..6cfe8f8 100644 --- a/examples/test_run.jl +++ b/examples/test_run.jl @@ -1,80 +1,48 @@ using JLD2, Random, Statistics, Subzero, BenchmarkTools, StructArrays import LibGEOS as LG -frac_stress = [-29955.396 -3428.008; -3428.008 -1942.0464] -frac_deform_floe = Floe( - [[ - [-50548.186, -49995.968], - [-50550.745, -37790.078], - [-20856.010, -32518.566], - [-20929.577, -49989.757], - [-50548.186, -49995.968], - ]], - 0.25, - 0.0, - u = 0.1, - v = -0.2, - ξ = 0.05, +Δt = 10 +max_overlap = 0.75 +Lx = 1e5 +Ly = Lx +hmean = 0.25 +Δh = 0.0 +grid = RegRectilinearGrid( + (-Lx, Lx), + (-Ly, Ly), + 1e4, + 1e4, ) -frac_floe = deepcopy(frac_deform_floe) # Without interactions, won't deform -no_frac_floe = Floe( # This floe is colliding with frac_deform_floe - [[ - [1467.795, -25319.563], - [1664.270, -25640.216], - [-1105.179, -33458.936], - [-17529.019, -50035.583], - [-21193.828, -50088.777], - [-21370.170, -32618.322], - [-21247.656, -31077.536], - [-12818.593, -27031.048], - [1467.795, -25319.563], - ]], - 0.25, - 0.0, +nboundary = PeriodicBoundary(North, grid) +sboundary = PeriodicBoundary(South, grid) +eboundary = CollisionBoundary(East, grid) +wboundary = OpenBoundary(West, grid) +topo_elem = TopographyElement( + [[[1e4, 0.0], [0.0, 1e4], [1e4, 2e4], [2e4, 1e4], [1e4, 0.0]]], +) +domain = Domain( + nboundary, + sboundary, + eboundary, + wboundary, + StructArray([topo_elem]), ) -no_frac_small = Floe( # This floe is too small to fracture or deform +# Diagonal floe barely overlaping with eastern collision boundary +efloe_small = Floe( [[ - [1e3, 1e3], - [1e3, 1.5e3], - [1.5e3, 1.5e3], - [1.5e3, 1e3], - [1e3, 1e3], + [9.5e4, 0.0], + [9e4, 0.5e4], + [10e4, 2.5e4], + [10.05e4, 2e4], + [9.5e4, 0.0], ]], - 0.25, - 0.0, -) -frac_deform_floe.stress = frac_stress -frac_deform_floe.interactions = collect([ - 3, - -279441968.984, - -54223517.438, - -21091.0918258529, - -40358.0042297616, - -148920620521.112, - 6795329.38154967, -]') -frac_deform_floe.p_dudt = 0.11 -frac_floe.stress = frac_stress -no_frac_small.stress = frac_stress - -floes = StructArray([ - frac_deform_floe, frac_floe, no_frac_floe, no_frac_small -]) -floes.id .= collect(1:4) -frac_settings = FractureSettings( - fractures_on = true, - criteria = HiblerYieldCurve(floes), - Δt = 75, - deform_on = true, + hmean, + Δh, ) -new_floes = Subzero.split_floe( - floes[1], - Xoshiro(3), - frac_settings, - CouplingSettings(), - Constants(), - 10, -) +efloe_small.u = 0.5 +efloe_small.v = 0.25 +consts = Constants() +Subzero.floe_domain_interaction!(efloe_small, domain, consts, Δt, max_overlap) # User Inputs const FT = Float64 diff --git a/src/floe_utils.jl b/src/floe_utils.jl index e72cf38..afb1972 100644 --- a/src/floe_utils.jl +++ b/src/floe_utils.jl @@ -54,7 +54,6 @@ Output: """ find_poly_coords(poly::LG.Polygon) = LG.GeoInterface.coordinates(poly)::PolyVec{Float64} -# [LG.GeoInterface.coordinates(LG.exteriorRing(poly))]::PolyVec{Float64} """ get_polygons(geom) @@ -65,36 +64,110 @@ Inputs: Outputs: """ -get_polygons(geom) = Vector{LG.Polygon}() +function get_polygons(geom) + polys = + if geom isa LG.Polygon + LG.area(geom) > 0 ? [geom] : Vector{LG.Polygon}() + elseif geom isa LG.MultiPolygon + sub_geoms = LG.getGeometries(geom) + for i in reverse(eachindex(sub_geoms)) + if LG.area(sub_geoms[i]) == 0 + deleteat!(sub_geoms, i) + end + end + sub_geoms + else + Vector{LG.Polygon}() + end + return polys::Vector{LG.Polygon} +end """ - get_polygons(poly::LG.Polygon) = [poly] + intersect_polys(p1, p2) -Return a polygon list with polygon element provided. +Intersect two geometries and return a list of polygons resulting. Inputs: - poly -Outputs: - + p1 + p2 +Output: + Vector of LibGEOS Polygons """ -get_polygons(poly::LG.Polygon) = [poly] +intersect_polys(p1, p2) = get_polygons( + LG.intersection( + p1, + p2, + )::LG.Geometry, +)::Vector{LG.Polygon} """ - get_polygons(multipoly::LG.MultiPolygon) + intersect_coords(c1, c2) -Returns a list of polygons that make up multipolygon input +Intersect geometries using their coordinates to return list of resulting +polygons. Inputs: - multipoly -Outputs: - sub_geoms + c1 + c2 +Output: + Vector of LibGEOS Polygons """ -function get_polygons(multipoly::LG.MultiPolygon) - sub_geoms = LG.getGeometries(multipoly) - for i in reverse(eachindex(sub_geoms)) - if LG.area(sub_geoms[i]) == 0 - deleteat!(sub_geoms, i) - end - end - return sub_geoms::Vector{LG.Polygon} +intersect_coords(c1, c2) = intersect_polys( + LG.Polygon(c1), + LG.Polygon(c2), +)::Vector{LG.Polygon} + +""" + deepcopy_floe(floe::LazyRow{Floe{FT}}) + +Deepcopy of a floe by creating a new floe and copying all fields. +Inputs: + floe +Outputs: + New floe with floes that are equal in value. Any vector fields are copies so + they share values, but not referance. +""" +function deepcopy_floe(floe::LazyRow{Floe{FT}}) where {FT} + f = Floe{FT}( + centroid = copy(floe.centroid), + coords = translate(floe.coords, 0, 0), + height = floe.height, + area = floe.area, + mass = floe.mass, + rmax = floe.rmax, + moment = floe.moment, + angles = copy(floe.angles), + x_subfloe_points = copy(floe.x_subfloe_points), + y_subfloe_points = copy(floe.y_subfloe_points), + α = floe.α, + u = floe.u, + v = floe.v, + ξ = floe.ξ, + status = Status(floe.status.tag, copy(floe.status.fuse_idx)), + id = floe.id, + ghost_id = floe.ghost_id, + parent_ids = copy(floe.parent_ids), + ghosts = copy(floe.ghosts), + fxOA= floe.fxOA, + fyOA = floe.fyOA, + trqOA = floe.trqOA, + hflx_factor = floe.hflx_factor, + overarea = floe.overarea, + collision_force = copy(floe.collision_force), + collision_trq = floe.collision_trq, + stress = copy(floe.stress), + stress_history = StressCircularBuffer{FT}( + capacity(floe.stress_history.cb), + ), + strain = copy(floe.strain), + p_dxdt = floe.p_dxdt, + p_dydt = floe.p_dydt, + p_dudt = floe.p_dudt, + p_dvdt = floe.p_dvdt, + p_dξdt = floe.p_dξdt, + p_dαdt = floe.p_dαdt, + ) + f.stress_history.total .= copy(floe.stress_history.total) + append!(f.stress_history.cb, copy(floe.stress_history.cb)) + return f end """ @@ -133,12 +206,14 @@ Output: Updates given coords """ function translate(coords::PolyVec{FT}, Δx, Δy) where {FT<:AbstractFloat} - new_coords = deepcopy(coords) - translate!(new_coords, Δx, Δy) + new_coords = [[Vector{Float64}(undef, 2) for _ in eachindex(coords[1])]] + for i in eachindex(coords[1]) + new_coords[1][i][1] = coords[1][i][1] + Δx + new_coords[1][i][2] = coords[1][i][2] + Δy + end return new_coords end - """ translate!(coords, Δx, Δy) @@ -708,56 +783,47 @@ Inputs: l1 line/polygon coordinates l2 line/polygon coordinates Outputs: - N intersection points in a Nx2 matrix where column 1 - is the x-coordinates and column 2 is the y-coordinates and each row is - an intersection point. - -Note - Translated into Julia from the following program: - NS (2022). Curve intersections - (https://www.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections), - MATLAB Central File Exchange. Retrieved November 2, 2022. - Only translated for the case where l1 and l2 are distinct. -""" -function intersect_lines(l1, l2) - x1, y1 = separate_xy(l1) - x2, y2 = separate_xy(l2) - x2t = x2' - y2t = y2' - Δx1 = diff(x1) - Δx2 = diff(x2t, dims = 2) - Δy1 = diff(y1) - Δy2 = diff(y2t, dims = 2) - - # Determine 'signed distances' - S1 = Δx1 .* @view(y1[1:end-1]) .- Δy1 .* @view(x1[1:end-1]) - s1 = (Δx1 .* y2t .- Δy1 .*x2t) # Needed for S1 calculation - C1 = (@view(s1[:, 1:end-1]) .- S1) .* (@view(s1[:, 2:end]) .- S1) .<= 0 - - S2 = Δx2 .* @view(y2t[:, 1:end-1]) .- Δy2 .* @view(x2t[:, 1:end-1]) - s2 = (y1 .* Δx2 .- x1 .* Δy2)' # Needed for S2 calculation - C2 = ((@view(s2[:, 1:end-1]) .- S2') .* (@view(s2[:, 2:end]) .- S2') .<= 0)' - - # Obtain the segments where an intersection is expected - idx = findall(C1 .& C2) - - P = if isempty(idx) - zeros(eltype(x1), 2,0) - else - # Transpose and prepare for output - i = getindex.(idx, 1) - j = getindex.(idx, 2)[:, :] - Δx2t = Δx2' - Δy2t = Δy2' - S2t = S2' - L = Δy2t[j] .* Δx1[i] - Δy1[i] .* Δx2t[j] - i = i[:, :][L .!= 0] - j = j[L .!= 0] - L = L[L .!= 0] - # Solve system of eqs to get the common points - unique(hcat((Δx2t[j] .* S1[i] - Δx1[i] .* S2t[j]) ./ L, - (Δy2t[j] .* S1[i] - Δy1[i] .* S2t[j]) ./ L), dims = 1) + Set of points that are at the intersection of the + two line segments. +""" +function intersect_lines(l1::PolyVec{FT}, l2) where {FT} + points = Set{Tuple{FT, FT}}() + for i in 1:(length(l1[1]) - 1) + # First line represented as a1x + b1y = c1 + x11 = l1[1][i][1] + x12 = l1[1][i+1][1] + y11 = l1[1][i][2] + y12 = l1[1][i+1][2] + a1 = y12 - y11 + b1 = x11 - x12 + c1 = a1 * x11 + b1 * y11 + for j in 1:(length(l2[1]) - 1) + # Second line represented as a2x + b2y = c2 + x21 = l2[1][j][1] + x22 = l2[1][j+1][1] + y21 = l2[1][j][2] + y22 = l2[1][j+1][2] + a2 = y22 - y21 + b2 = x21 - x22 + c2 = a2 * x21 + b2 * y21 + + determinant = a1 * b2 - a2 * b1 + # Find place there two lines cross + # Note that lines extend beyond given line segments + if determinant != 0 + x = (b2*c1 - b1*c2)/determinant + y = (a1*c2 - a2*c1)/determinant + # Make sure intersection is on given line segments + if min(x11, x12) <= x <= max(x11, x12) && + min(x21, x22) <= x <= max(x21, x22) && + min(y11, y12) <= y <= max(y11, y12) && + min(y21, y22) <= y <= max(y21, y22) + push!(points, (x, y)) + end + end + end end - return P + return points end """ diff --git a/src/output.jl b/src/output.jl index 637a202..b923df2 100644 --- a/src/output.jl +++ b/src/output.jl @@ -482,14 +482,23 @@ Output: function write_data!(sim, tstep) # write initial state on first timestep if tstep == 0 - write_init_state_data!(sim, tstep) + write_init_state_data!(sim) end # Write checkpoint data write_checkpoint_data!(sim, tstep) # Write floe data - write_floe_data!(sim, tstep) + write_floe_data!( + sim.writers.floewriters, + sim.model.floes, + tstep, + ) # Write grid data - write_grid_data!(sim, tstep) + write_grid_data!( + sim.writers.gridwriters, + sim.model.floes, + sim.model.domain.topography, + tstep, + ) return end @@ -503,7 +512,7 @@ Inputs: Outputs: Saves simulation state to file. """ -function write_init_state_data!(sim, tstep) +function write_init_state_data!(sim) for filepath in sim.writers.initialwriters.filepath jldopen(filepath, "a") do file file["sim"] = sim @@ -543,45 +552,48 @@ end Writes desired FloeOutputWriter data to JLD2 file. Inputs: - sim simulation to run + writers list of floe writers + floes list of floes tstep simulation timestep Output: Writes desired fields writer.outputs to JLD2 file with name writer.fn for current timestep, which will be the group in the JLD2 file. """ -function write_floe_data!(sim, tstep) - for w in sim.writers.floewriters[ - mod.(tstep, sim.writers.floewriters.Δtout) .== 0 +function write_floe_data!(writers, floes, tstep) + for w in writers[ + mod.(tstep, writers.Δtout) .== 0 ] - jldopen(w.filepath, "a+") do file - for output in w.outputs - file[string(output, "/", tstep)] = - StructArrays.component(sim.model.floes, output) - end + file = jldopen(w.filepath, "a+") + for output in w.outputs + vals = StructArrays.component(floes, output) + write_field(file, string(output, "/", tstep), vals) end + close(file) end # once parents are recorded, they should be deleted - empty!.(sim.model.floes.parent_ids) + empty!.(floes.parent_ids) return end +write_field(file, group, vals) = write(file, group, vals) """ write_grid_data!(sim, tstep) Writes desired GridOutputWriter data to NetCDF file. Inputs: - sim simulation to run + writers list of grid writers + floes list of floes + topography list of topography elements tstep simulation timestep Output: Writes desired fields writer.outputs to file with name writer.fn for current timestep. """ -function write_grid_data!(sim, tstep) - live_floes = filter(f -> f.status.tag == active, sim.model.floes) - w_idx = mod.(tstep, sim.writers.gridwriters.Δtout) .== 0 - if !isempty(live_floes) && sum(w_idx) != 0 - for w in sim.writers.gridwriters[w_idx] - calc_eulerian_data!(live_floes, sim.model.domain.topography, w) +function write_grid_data!(writers, floes, topography, tstep) + w_idx = mod.(tstep, writers.Δtout) .== 0 + if !isempty(floes) && sum(w_idx) != 0 + @views for w in writers[w_idx] + calc_eulerian_data!(floes, topography, w) istep = div(tstep, w.Δtout) + 1 # Julia indicies start at 1 # Open file and write data from grid writer @@ -839,68 +851,69 @@ function calc_eulerian_data!(floes, topography, writer) pic_polys = [ LG.intersection( cell_poly, - LG.Polygon(floes[idx].coords), + LG.Polygon(floes.coords[idx]), ) for idx in floeidx] pic_area = [LG.area(poly) for poly in pic_polys] floeidx = floeidx[pic_area .> 0] pic_area = pic_area[pic_area .> 0] fic = floes[floeidx] - fic_area = fic.area + fic_area = floes.area[floeidx] + fic_mass = floes.mass[floeidx] area_ratios = pic_area ./ fic_area area_tot = sum(pic_area) - mass_tot = sum(fic.mass .* area_ratios) + mass_tot = sum(floes.mass[floeidx] .* area_ratios) if mass_tot > 0 # mass and area ratios - ma_ratios = area_ratios .* (fic.mass ./ mass_tot) + ma_ratios = area_ratios .* (floes.mass[floeidx] ./ mass_tot) outputs = writer.outputs for k in eachindex(outputs) data = if outputs[k] == :u_grid - sum(fic.u .* ma_ratios) + sum(floes.u[floeidx] .* ma_ratios) elseif outputs[k] == :v_grid - sum(fic.v .* ma_ratios) + sum(floes.v[floeidx] .* ma_ratios) elseif outputs[k] == :dudt_grid - sum(fic.p_dudt .* ma_ratios) + sum(floes.p_dudt[floeidx] .* ma_ratios) elseif outputs[k] == :dvdt_grid - sum(fic.p_dvdt .* ma_ratios) + sum(floes.p_dvdt[floeidx] .* ma_ratios) elseif outputs[k] == :si_frac_grid area_tot/LG.area(cell_poly) elseif outputs[k] == :overarea_grid - sum(fic.overarea)/length(fic) + sum(floes.overarea[floeidx])/length(floeidx) elseif outputs[k] == :mass_grid mass_tot elseif outputs[k] == :area_grid area_tot elseif outputs[k] == :height_grid - sum(fic.height .* ma_ratios) + sum(floes.height[floeidx] .* ma_ratios) elseif outputs[k] == :stress_xx_grid - sum([s[1, 1] for s in fic.stress] .* ma_ratios) + sum([s[1, 1] for s in floes.stress[floeidx]] .* ma_ratios) elseif outputs[k] == :stress_yx_grid - sum([s[1, 2] for s in fic.stress] .* ma_ratios) + sum([s[1, 2] for s in floes.stress[floeidx]] .* ma_ratios) elseif outputs[k] == :stress_xy_grid - sum([s[2, 1] for s in fic.stress] .* ma_ratios) + sum([s[2, 1] for s in floes.stress[floeidx]] .* ma_ratios) elseif outputs[k] == :stress_yy_grid - sum([s[2, 2] for s in fic.stress] .* ma_ratios) + sum([s[2, 2] for s in floes.stress[floeidx]] .* ma_ratios) elseif outputs[k] == :stress_eig_grid - xx = sum([s[1, 1] for s in fic.stress] .* ma_ratios) - yx = sum([s[1, 2] for s in fic.stress] .* ma_ratios) - xy = sum([s[2, 1] for s in fic.stress] .* ma_ratios) - yy = sum([s[2, 2] for s in fic.stress] .* ma_ratios) + xx = sum([s[1, 1] for s in floes.stress[floeidx]] .* ma_ratios) + yx = sum([s[1, 2] for s in floes.stress[floeidx]] .* ma_ratios) + xy = sum([s[2, 1] for s in floes.stress[floeidx]] .* ma_ratios) + yy = sum([s[2, 2] for s in floes.stress[floeidx]] .* ma_ratios) stress = maximum(eigvals([xx yx; xy yy])) if abs(stress) > 1e8 stress = 0.0 end stress elseif outputs[k] == :strain_ux_grid - sum([s[1, 1] for s in fic.strain] .* ma_ratios) + sum([s[1, 1] for s in floes[floeidx].strain] .* ma_ratios) elseif outputs[k] == :strain_vx_grid - sum([s[1, 2] for s in fic.strain] .* ma_ratios) + sum([s[1, 2] for s in floes[floeidx].strain] .* ma_ratios) elseif outputs[k] == :strain_uy_grid - sum([s[2, 1] for s in fic.strain] .* ma_ratios) + sum([s[2, 1] for s in floes[floeidx].strain] .* ma_ratios) elseif outputs[k] == :strain_vy_grid - sum([s[2, 2] for s in fic.strain] .* ma_ratios) + sum([s[2, 2] for s in floes[floeidx].strain] .* ma_ratios) end writer.data[j, i, k] = data end diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl index dea17b2..ed1eee2 100644 --- a/src/physical_processes/collisions.jl +++ b/src/physical_processes/collisions.jl @@ -25,7 +25,6 @@ Inputs: and 2 where the first column is the x-coordinates and the second column is the y-coordinates force_factor Spring constant equivalent for collisions - t Float type model is running on (Float64 or Float32) Outputs: normal force of collision Δl mean length of distance between intersection points @@ -34,32 +33,41 @@ function calc_normal_force( c1, c2, region, - area::FT, + area, ipoints, force_factor::FT, ) where {FT<:AbstractFloat} force_dir = zeros(FT, 2) coords = find_poly_coords(region) - n_ipoints = size(ipoints, 1) + n_ipoints = length(ipoints) # Identify which region coordinates are the intersection points (ipoints) - verts = zeros(Int64, n_ipoints) - dists = zeros(FT, n_ipoints) - for i in 1:n_ipoints - p_i = repeat(ipoints[i, :], length(coords)) - dists[i], verts[i] = findmin([sum((c .- p_i).^2) for c in coords[1]]) + #verts = zeros(Int64, n_ipoints) + #dists = zeros(FT, n_ipoints) + p = Vector{Vector{FT}}() + @views for ip in ipoints + min_dist = FT(Inf) + min_vert = 1 + for i in eachindex(coords[1]) + dist = sqrt((coords[1][i][1] - ip[1])^2 + (coords[1][i][2] - ip[2])^2) + if dist < min_dist + min_dist = dist + min_vert = i + end + end + if min_dist < 1 + push!(p, coords[1][min_vert]) + end end - dists = sqrt.(dists) - p = coords[1][verts[findall(d -> d<1, dists)]] # maybe do rmax/1000 m = length(p) # Calculate force direction Δl = FT(0) if m == 2 # Two intersection points - Δx = p[2][1] - p[1][1] - Δy = p[2][2] - p[1][2] + Δx = FT(p[2][1] - p[1][1]) + Δy = FT(p[2][2] - p[1][2]) Δl = sqrt(Δx^2 + Δy^2) if Δl > 0.1 # should match our scale - force_dir = [-Δy/Δl; Δx/Δl] + force_dir .= [-Δy/Δl; Δx/Δl] end elseif m != 0 # Unusual number of intersection points x, y = separate_xy(coords) @@ -80,23 +88,23 @@ function calc_normal_force( Δl = mean(mag[on_idx]) if Δl > 0.1 Fn_tot = sum(Fn[on_idx, :], dims = 1) - force_dir = Fn_tot'/sqrt(Fn_tot*Fn_tot') + force_dir .= Fn_tot'/sqrt(Fn_tot*Fn_tot')[1] end end end # Check if direction of the force desceases overlap, else negate direction if Δl > 0.1 - c1new = [[c .+ vec(force_dir) for c in c1[1]]] + c1new = translate(c1, force_dir[1], force_dir[2]) # Floe/boudary intersection after being moved in force direction - new_inter_floe = LG.intersection(LG.Polygon(c1new), LG.Polygon(c2)) + new_regions_list = intersect_coords(c1new, c2) # See if the area of overlap has increased in corresponding region - for new_region in LG.getGeometries(new_inter_floe) + for new_region in new_regions_list if LG.intersects(new_region, region) && LG.area(new_region)/area > 1 - force_dir *= -1 + force_dir .*= -1 end end end - return (force_dir * area * force_factor)', Δl + return force_dir * area * force_factor, Δl end """ @@ -131,25 +139,29 @@ function calc_elastic_forces( c1, c2, regions, - region_areas::Vector{FT}, + region_areas, force_factor::FT, ) where {FT<:AbstractFloat} ipoints = intersect_lines(c1, c2) # Intersection points ncontact = 0 - overlap = FT(0) - if !isempty(ipoints) && size(ipoints,2) >= 2 + if !isempty(ipoints) && length(ipoints) >= 2 # Find overlapping regions greater than minumum area n1 = length(c1[1]) - 1 n2 = length(c2[1]) - 1 min_area = min(n1, n2) * 100 / 1.75 - regions = regions[region_areas .> min_area] - region_areas = region_areas[region_areas .> min_area] - overlap = region_areas - ncontact = length(regions) + for i in reverse(eachindex(region_areas)) + if region_areas[i] < min_area + deleteat!(region_areas, i) + deleteat!(regions, i) + else + ncontact += 1 + end + end end # Calculate forces for each remaining region force = zeros(FT, ncontact, 2) fpoint = zeros(FT, ncontact, 2) + overlap = ncontact > 0 ? region_areas : zeros(FT, ncontact) Δl_lst = zeros(FT, ncontact) for k in 1:ncontact if region_areas[k] != 0 @@ -164,7 +176,8 @@ function calc_elastic_forces( ipoints, force_factor ) - force[k, :] = normal_force + force[k, 1] = normal_force[1] + force[k, 2] = normal_force[2] Δl_lst[k] = Δl end end @@ -172,46 +185,113 @@ function calc_elastic_forces( end """ - calc_friction_forces(v1, v2, fpoint, nforce, Δl, consts) + get_velocity( + floe, + x, + y, + ) -Calculate frictional force for collision between two floes. +Get velocity of a point, assumed to be on given floe. +Inputs: + floe floe + x x-coordinate of point to find velocity at + y y-coordinate of point to find velocity at +Outputs: + u u velocity at point (x, y) assuming it is on given floe + v v velocity at point (x, y) assuming it is on given floe +""" +function get_velocity( + floe::Union{LazyRow{Floe{FT}}, Floe{FT}}, + x::FT, + y::FT, +) where {FT} + u = floe.u + floe.ξ * (x - floe.centroid[1]) + v = floe.v + floe.ξ * (y - floe.centroid[2]) + return u, v +end + +""" + get_velocity( + element::AbstractDomainElement{FT}, + x, + y, + ) + +Get velocity, which is 0m/s by default, of a point on topography element or +boundary. +""" +get_velocity( + element::AbstractDomainElement{FT}, + x, + y, +) where {FT} = + FT(0), FT(0) + +""" + calc_friction_forces( + ifloe, + jfloe, + fpoints, + normal::Matrix{FT}, + Δl, + consts, + Δt, + ) +Calculate frictional force for collision between two floes or a floe and a +domain element. Input: - v1 Matrix of floe speeds for first floe in collision - at each point a force is applied to the floe - row is [u v] and one - row per each N point - v2 vector of floe speeds for second floe or boundary - in collision - see v1 for form - fpoint x,y-coordinates of the point the force is + ifloe first floe in collsion + jfloe either second floe or topography + element/boundary element + fpoints x,y-coordinates of the point the force is applied on floe overlap region normal x,y normal force applied on fpoint on floe overlap region - Δl mean length of distance between intersection points + Δl mean length of distance between intersection points consts model constants needed for calculations + Δt simulation's timestep Outputs: - frictional/tangential force of the collision + force frictional/tangential force of the collision in + x and y (each row) for each collision (each column) """ function calc_friction_forces( - v1, - v2, + ifloe, + jfloe, # could be a boundary or a topography + fpoints, normal::Matrix{FT}, Δl, consts, Δt, ) where {FT} - force = zeros(FT, size(v1, 1), 2) - G = consts.E/(2*(1+consts.ν)) # Shear modulus - # Difference in velocities between floes in x and y direction - vdiff = v1 .- v2 - # Friction forces for each vector - for i in axes(vdiff, 1) - v = vdiff[i, :] - n = normal[i, :] - force_dir = maximum(abs.(v)) == 0 ? zeros(FT, 2) : v/norm(v) - friction = G * Δl[i] * Δt * norm(n) * -dot(force_dir, v) * force_dir - if norm(friction) > consts.μ*norm(n) - friction = -consts.μ*norm(n)*force_dir + force = zeros(FT, size(fpoints, 1), 2) + G = consts.E/(2*(1+consts.ν)) + for i in axes(fpoints, 1) + px = fpoints[i, 1] + py = fpoints[i, 2] + nnorm = sqrt(normal[i, 1]^2 + normal[i, 2]^2) + iu, iv = get_velocity(ifloe, px, py) + ju, jv = get_velocity(jfloe, px, py) + + udiff = iu - ju + vdiff = iv - jv + + vnorm = sqrt(udiff^2 + vdiff^2) + xdir = FT(0) + ydir = FT(0) + if udiff != 0 || vdiff != 0 + xdir = udiff/vnorm + ydir = vdiff/vnorm + end + dot_dir = xdir * udiff + ydir * vdiff + xfriction = G * Δl[i] * Δt * nnorm * xdir * -dot_dir + yfriction = G * Δl[i] * Δt * nnorm * ydir * -dot_dir + norm_fric = sqrt(xfriction^2 + yfriction^2) + if norm_fric > consts.μ * nnorm + xfriction = -consts.μ * nnorm * xdir + yfriction = -consts.μ * nnorm * ydir end - force[i, :] = friction + force[i, 1] = xfriction + force[i, 2] = yfriction end return force end @@ -262,19 +342,13 @@ function floe_floe_interaction!( Δt, max_overlap::FT, ) where {FT<:AbstractFloat} - ifloe_poly = LG.Polygon(ifloe.coords) - jfloe_poly = LG.Polygon(jfloe.coords) - inter_floe = LG.intersection(ifloe_poly, jfloe_poly) - inter_regions = LG.Polygon[] - region_areas = FT[] + inter_regions = intersect_coords(ifloe.coords, jfloe.coords) + region_areas = Vector{FT}(undef, length(inter_regions)) total_area = FT(0) - for geom in LG.getGeometries(inter_floe) - a = LG.area(geom)::FT - if a > 0 - push!(inter_regions, geom) - push!(region_areas, a) - total_area += a - end + for i in eachindex(inter_regions) + a = FT(LG.area(inter_regions[i])) + region_areas[i] = a + total_area += a end if total_area > 0 # Floes overlap too much - remove floe or transfer floe mass @@ -285,9 +359,9 @@ function floe_floe_interaction!( else # Constant needed to calculate force ih = ifloe.height - ir = sqrt(ifloe.area) + ir = FT(sqrt(ifloe.area)) jh = jfloe.height - jr = sqrt(jfloe.area) + jr = FT(sqrt(jfloe.area)) force_factor = if ir>1e5 || jr>1e5 consts.E*min(ih, jh)/min(ir, jr) else @@ -305,23 +379,25 @@ function floe_floe_interaction!( velocities at force points =# np = size(fpoints, 1) if np > 0 - iv = repeat([ifloe.u ifloe.v], outer = np) .+ - ifloe.ξ*(fpoints .- repeat(ifloe.centroid', outer = np)) - jv = repeat([jfloe.u jfloe.v], outer = np) .+ - jfloe.ξ*(fpoints .- repeat(jfloe.centroid', outer = np)) friction_forces = calc_friction_forces( - iv, - jv, + ifloe, + jfloe, + fpoints, normal_forces, Δl, consts, - Δt + Δt, ) # Calculate total forces and update ifloe's interactions forces = normal_forces .+ friction_forces if sum(abs.(forces)) != 0 - ifloe.interactions = [ifloe.interactions; - fill(j, np) forces fpoints zeros(np) overlaps] + new_interactions = [ + fill(j, np) forces fpoints zeros(FT, np) overlaps + ] + ifloe.interactions = vcat( + ifloe.interactions, + new_interactions + ) ifloe.overarea += sum(overlaps) end end @@ -533,23 +609,17 @@ function floe_domain_element_interaction!( Δt, max_overlap::FT, ) where {FT} - floe_poly = LG.Polygon(floe.coords) - bounds_poly = LG.Polygon(element.coords) - # Check if the floe and element actually overlap - inter_floe = LG.intersection(floe_poly, bounds_poly) - inter_regions = LG.Polygon[] - region_areas = FT[] + inter_regions = intersect_coords(floe.coords, element.coords) + region_areas = Vector{FT}(undef, length(inter_regions)) max_area = FT(0) - for geom in LG.getGeometries(inter_floe) - a = LG.area(geom)::FT - if a > 0 - push!(inter_regions, geom) - push!(region_areas, a) - if a > max_area - max_area = a - end + for i in eachindex(inter_regions) + a = FT(LG.area(inter_regions[i])) + region_areas[i] = a + if a > max_area + max_area = a end end + if max_area > 0 # Regions overlap too much if maximum(region_areas)/floe.area > max_overlap @@ -568,23 +638,23 @@ function floe_domain_element_interaction!( normal_direction_correct!(normal_forces, fpoints, element) # Calculate frictional forces at each force point np = size(fpoints, 1) - vfloe = repeat([floe.u floe.v], outer = np) .+ - floe.ξ*(fpoints .- repeat(floe.centroid', outer = np)) - vbound = repeat(zeros(FT, 1, 2), outer = np) - friction_forces = calc_friction_forces( - vfloe, - vbound, - normal_forces, - Δl, - consts, - Δt, - ) - # Calculate total forces and update ifloe's interactions - forces = normal_forces .+ friction_forces - if sum(abs.(forces)) != 0 - floe.interactions = [floe.interactions; - fill(Inf, np) forces fpoints zeros(np) overlaps] - floe.overarea += sum(overlaps) + if np > 0 + friction_forces = calc_friction_forces( + floe, + element, + fpoints, + normal_forces, + Δl, + consts, + Δt, + ) + # Calculate total forces and update ifloe's interactions + forces = normal_forces .+ friction_forces + if sum(abs.(forces)) != 0 + floe.interactions = [floe.interactions; + fill(Inf, np) forces fpoints zeros(np) overlaps] + floe.overarea += sum(overlaps) + end end end end @@ -752,7 +822,9 @@ Inputs: Outputs: None. Floe's interactions field is updated with calculated torque. """ -function calc_torque!(floe::Union{LazyRow{<:Floe{FT}}, Floe{FT}}) where {FT<:AbstractFloat} +function calc_torque!( + floe::Union{LazyRow{<:Floe{FT}}, Floe{FT}}, +) where {FT<:AbstractFloat} inters = floe.interactions if !isempty(inters) dir = [inters[:, xpoint] .- floe.centroid[1] inters[:, ypoint] .- floe.centroid[2] zeros(FT, size(inters, 1))] @@ -764,18 +836,33 @@ function calc_torque!(floe::Union{LazyRow{<:Floe{FT}}, Floe{FT}}) where {FT<:Abs floe.interactions[i, torque] = itorque[3] end end + return end """ - + potential_interaction( + centroid1, + centroid2, + rmax1, + rmax2, + ) +Determine if two floes could potentially interact using the two centroid and two +radii to form a bounding circle. +Inputs: + centroid1 first floe's centroid [x, y] + centroid2 second floe's centroid [x, y] + rmax1 first floe's maximum radius + rmax2 second floe's maximum radius +Outputs: + true if floes could potentially interact, false otherwise """ potential_interaction( centroid1, centroid2, rmax1, rmax2, -) = (sum((centroid1 .- centroid2).^2) < (rmax1 + rmax2)^2) - +) = ((centroid1[1] - centroid2[1])^2 + (centroid1[2] - centroid2[2])^2) < + (rmax1 + rmax2)^2 """ timestep_collisions!( @@ -810,7 +897,7 @@ function timestep_collisions!( ) where {FT<:AbstractFloat} collide_pairs = Dict{Tuple{Int, Int}, Tuple{Int, Int}}() # floe-floe collisions for floes i and j where i i jidx = Int(j) - floes.interactions[jidx] = [floes.interactions[jidx]; i -ij_inters[inter_idx, xforce] -ij_inters[inter_idx, yforce] #= - =# ij_inters[inter_idx, xpoint] ij_inters[inter_idx, ypoint] FT(0.0) ij_inters[inter_idx, overlap]] + floes.interactions[jidx] = vcat( + floes.interactions[jidx], + ij_inters[inter_idx:inter_idx, :] + ) + floes.interactions[jidx][end, floeidx] = i + floes.interactions[jidx][end, xforce] *= -1 + floes.interactions[jidx][end, yforce] *= -1 + # floes.interactions[jidx] = [floes.interactions[jidx]; i -ij_inters[inter_idx, xforce] -ij_inters[inter_idx, yforce] #= + # =# ij_inters[inter_idx, xpoint] ij_inters[inter_idx, ypoint] FT(0.0) ij_inters[inter_idx, overlap]] floes.overarea[jidx] += ij_inters[inter_idx, overlap] end end @@ -908,31 +1002,48 @@ end If the given element intersects with the boundary, add ghosts of the element and any of its existing ghosts. Inputs: - element given element - ghosts current ghosts of element - boundary boundary to translate element through - trans_vec 1x2 matrix of form [x y] to translate element through the boundary + floes model's list of floes + elem_idx floe of interest's index within the floe list + boundary boundary to translate element through + trans_vec 1x2 matrix of form [x y] to translate element + through the boundary Outputs: - New ghosts created by the given element, or its current ghosts, passing through the given boundary. + Nothing. New ghosts created by the given element, or its current ghosts, + are added to the floe list """ -function ghosts_on_bounds(element, ghosts, boundary, trans_vec) - new_ghosts = StructArray(Vector{typeof(element)}()) - if LG.area(LG.intersection(LG.Polygon(element.coords), LG.Polygon(boundary.coords))) > 0 +function ghosts_on_bounds!( + floes, + elem_idx, + boundary, + trans_vec, +) + nfloes = length(floes) + nghosts = 1 + if !isempty(intersect_coords(floes.coords[elem_idx], boundary.coords)) # ghosts of existing ghosts and original element - append!(new_ghosts, deepcopy.(ghosts)) - push!(new_ghosts, deepcopy(element)) - for i in eachindex(new_ghosts) - translate!(new_ghosts.coords[i], trans_vec[1], trans_vec[2]) - new_ghosts.centroid[i] .+= trans_vec + for i in floes.ghosts[elem_idx] + push!( + floes, + deepcopy_floe(LazyRow(floes, i)), + ) + nghosts += 1 + end + push!( + floes, + deepcopy_floe(LazyRow(floes, elem_idx)), + ) + for i in (nfloes + 1):(nfloes + nghosts) + translate!(floes.coords[i], trans_vec[1], trans_vec[2]) + floes.centroid[i] .+= trans_vec end end - return new_ghosts + return end """ find_ghosts( - elem, - current_ghosts, + floes, + elem_idx, ebound::PeriodicBoundary, wbound::PeriodicBoundary, ) @@ -942,43 +1053,61 @@ periodic boundary. If element's centroid isn't within the domain in the east/west direction, swap it with its ghost since the ghost's centroid must then be within the domain. Inputs: - elem given element - current_ghosts current ghosts of element - eboundary domain's eastern boundary - wboundary domain's western boundary + floes model's list of floes + elem_idx floe of interest's index within the floe list + eboundary domain's eastern boundary + wboundary domain's western boundary Outputs: - Return "primary" element, which has its centroid within the domain in the - east/west direction, and all of its ghosts in the east/west direction, - including ghosts of previously existing ghosts. + None. Ghosts added to the floe list. Primary floe always has centroid within + the domain, else it is swapped with one of its ghost's which has a centroid + within the domain. """ function find_ghosts( - elem, - current_ghosts, + floes, + elem_idx, ebound::PeriodicBoundary{East, FT}, wbound::PeriodicBoundary{West, FT}, ) where {FT <: AbstractFloat} Lx = ebound.val - wbound.val - new_ghosts = - # passing through western boundary - if (elem.centroid[1] - elem.rmax < wbound.val) - ghosts_on_bounds(elem, current_ghosts, wbound, [Lx, FT(0)]) - # passing through eastern boundary - elseif (elem.centroid[1] + elem.rmax > ebound.val) - ghosts_on_bounds(elem, current_ghosts, ebound, [-Lx, FT(0)]) - else - StructArray(Vector{typeof(elem)}()) - end + nfloes = length(floes) + # passing through western boundary + if (floes.centroid[elem_idx][1] - floes.rmax[elem_idx] < wbound.val) + ghosts_on_bounds!( + floes, + elem_idx, + wbound, + [Lx, FT(0)], + ) + # passing through eastern boundary + elseif (floes.centroid[elem_idx][1] + floes.rmax[elem_idx] > ebound.val) + ghosts_on_bounds!( + floes, + elem_idx, + ebound, + [-Lx, FT(0)], + ) + end # if element centroid isn't in domain's east/west direction, swap with ghost - if !isempty(new_ghosts) && ((elem.centroid[1] < wbound.val) || (ebound.val < elem.centroid[1])) - elem, new_ghosts[end] = new_ghosts[end], elem + if length(floes) > nfloes + if floes.centroid[elem_idx][1] < wbound.val + translate!(floes.coords[elem_idx], Lx, FT(0)) + floes.centroid[elem_idx][1] += Lx + translate!(floes.coords[end], -Lx, FT(0)) + floes.centroid[end][1] -= Lx + elseif ebound.val < floes.centroid[elem_idx][1] + translate!(floes.coords[elem_idx], -Lx, FT(0)) + floes.centroid[elem_idx][1] -= Lx + translate!(floes.coords[end], Lx, FT(0)) + floes.centroid[end][1] += Lx + end end - return elem, new_ghosts + return end """ find_ghosts( - elem, - current_ghosts, + floes, + elem_idx, nbound::PeriodicBoundary{North, <:AbstractFloat}, sbound::PeriodicBoundary{South, <:AbstractFloat}, ) @@ -988,37 +1117,55 @@ southern periodic boundary. If element's centroid isn't within the domain in the north/south direction, swap it with its ghost since the ghost's centroid must then be within the domain. Inputs: - elem given element - current_ghosts current ghosts of element + floes model's list of floes + elem_idx floe of interest's index within the floe list nboundary domain's northern boundary sboundary domain's southern boundary Outputs: - Return "primary" element, which has its centroid within the domain in the - north/south direction, and all of its ghosts in the north/south direction, - including ghosts of previously existing ghosts. + None. Ghosts added to the floe list. Primary floe always has centroid within + the domain, else it is swapped with one of its ghost's which has a centroid + within the domain. """ function find_ghosts( - elem, - current_ghosts, + floes, + elem_idx, nbound::PeriodicBoundary{North, FT}, sbound::PeriodicBoundary{South, FT}, ) where {FT <: AbstractFloat} Ly = nbound.val - sbound.val - new_ghosts = + nfloes = length(floes) # passing through southern boundary - if (elem.centroid[2] - elem.rmax < sbound.val) - ghosts_on_bounds(elem, current_ghosts, sbound, [FT(0), Ly]) - # passing through northern boundary - elseif (elem.centroid[2] + elem.rmax > nbound.val) - ghosts_on_bounds(elem, current_ghosts, nbound, [FT(0), -Ly]) - else - StructArray(Vector{typeof(elem)}()) - end + if (floes.centroid[elem_idx][2] - floes.rmax[elem_idx] < sbound.val) + ghosts_on_bounds!( + floes, + elem_idx, + sbound, + [FT(0), Ly], + ) + # passing through northern boundary + elseif (floes.centroid[elem_idx][2] + floes.rmax[elem_idx] > nbound.val) + ghosts_on_bounds!( + floes, + elem_idx, + nbound, + [FT(0), -Ly], + ) + end # if element centroid isn't in domain's north/south direction, swap with ghost - if !isempty(new_ghosts) && ((elem.centroid[2] < sbound.val) || (nbound.val < elem.centroid[2])) - elem, new_ghosts[end] = new_ghosts[end], elem + if length(floes) > nfloes + if floes.centroid[elem_idx][2] < sbound.val + translate!(floes.coords[elem_idx], FT(0), Ly) + floes.centroid[elem_idx][2] += Ly + translate!(floes.coords[end], FT(0), -Ly) + floes.centroid[end][2] -= Ly + elseif nbound.val < floes.centroid[elem_idx][2] + translate!(floes.coords[elem_idx], FT(0), -Ly) + floes.centroid[elem_idx][2] -= Ly + translate!(floes.coords[end], FT(0), Ly) + floes.centroid[end][2] += Ly + end end - return elem, new_ghosts + return end """ @@ -1034,39 +1181,36 @@ Outputs: None. Ghosts of floes are added to floe list. """ function add_floe_ghosts!( - floes::StructArray{Floe{FT}}, + floes::FLT, max_boundary, min_boundary, -) where {FT <: AbstractFloat} +) where {FT <: AbstractFloat, FLT <: StructArray{<:Floe{FT}}} nfloes = length(floes) # uses initial length of floes so we can append to list for i in eachindex(floes) - f = floes[i] # the floe is active in the simulation and a parent floe - if f.status.tag == active && f.ghost_id == 0 - f, new_ghosts = find_ghosts( - f, - floes[f.ghosts], + if floes.status[i].tag == active && floes.ghost_id[i] == 0 + find_ghosts( + floes, + i, max_boundary, min_boundary, ) - if !isempty(new_ghosts) - nghosts = length(new_ghosts) - new_ghosts.ghost_id .= range(1, nghosts) .+ length(f.ghosts) + new_nfloes = length(floes) + if new_nfloes > nfloes + nghosts = new_nfloes - nfloes + # set ghost floe's ghost id + floes.ghost_id[(nfloes + 1):new_nfloes] .= range(1, nghosts) .+ length(floes.ghosts[i]) # remove ghost floes ghosts as these were added to parent - empty!.(new_ghosts.ghosts) - # add ghosts to floe list - append!(floes, new_ghosts) + empty!.(floes.ghosts[nfloes + 1:new_nfloes]) # index of ghosts floes saved in parent - append!(f.ghosts, (nfloes + 1):(nfloes + nghosts)) + append!(floes.ghosts[i], (nfloes + 1):(nfloes + nghosts)) nfloes += nghosts - floes[i] = f end end end return end - """ add_ghosts!( elems, diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index 216fc1f..d311c30 100644 --- a/src/physical_processes/coupling.jl +++ b/src/physical_processes/coupling.jl @@ -1580,7 +1580,7 @@ function calc_two_way_coupling!( ) where {FT} # Determine force from floe on each grid cell it is in cell_area = grid.Δx * grid.Δy - Threads.@threads for cartidx in CartesianIndices(ocean.scells) + for cartidx in CartesianIndices(ocean.scells) ocean.τx[cartidx] = FT(0) ocean.τy[cartidx] = FT(0) ocean.si_frac[cartidx] = FT(0) @@ -1603,10 +1603,13 @@ function calc_two_way_coupling!( floe_locations.Δy[i], ) floe_poly = LG.Polygon(floe_coords) - floe_area_in_cell = LG.area(LG.intersection( - cell_poly, - floe_poly, - ))::FT + floe_area_in_cell = FT(sum( + LG.area.(intersect_polys(cell_poly, floe_poly)) + )) + # floe_area_in_cell = FT(LG.area(LG.intersection( + # cell_poly, + # floe_poly, + # ))) if floe_area_in_cell > 0 # Add forces and area to ocean fields ocean.τx[cartidx] += (τocn.τx[i]/τocn.npoints[i]) * floe_area_in_cell diff --git a/src/physical_processes/update_floe.jl b/src/physical_processes/update_floe.jl index 05db082..160b16d 100644 --- a/src/physical_processes/update_floe.jl +++ b/src/physical_processes/update_floe.jl @@ -245,14 +245,11 @@ function calc_stress!(floe) inters = floe.interactions # Calculates timestep stress stress = fill(1/(2* floe.area * floe.height), 2, 2) - stress[1, 1] *= sum((inters[:, xpoint] .- xi) .* inters[:, xforce]) + - sum(inters[:, xforce] .* (inters[:, xpoint] .- xi)) + stress[1, 1] *= 2sum((inters[:, xpoint] .- xi) .* inters[:, xforce]) stress[1, 2] *= sum((inters[:, ypoint] .- yi) .* inters[:, xforce]) + sum(inters[:, yforce] .* (inters[:, xpoint] .- xi)) - stress[2, 1] *= sum((inters[:, xpoint] .- xi) .* inters[:, yforce]) + - sum(inters[:, xforce] .* (inters[:, ypoint] .- yi)) - stress[2, 2] *= sum((inters[:, ypoint] .- yi) .* inters[:, yforce]) + - sum(inters[:, yforce] .* (inters[:, ypoint] .- yi)) + stress[2, 1] = stress[1, 2] + stress[2, 2] *= 2sum((inters[:, ypoint] .- yi) .* inters[:, yforce]) # Add timestep stress to stress history push!(floe.stress_history, stress) # Average stress history to find floe's average stress diff --git a/src/simulation_components/domain_and_grid.jl b/src/simulation_components/domain_and_grid.jl index 993102b..4a8b495 100644 --- a/src/simulation_components/domain_and_grid.jl +++ b/src/simulation_components/domain_and_grid.jl @@ -585,25 +585,27 @@ struct Domain{ SB<:AbstractBoundary{South, FT}, EB<:AbstractBoundary{East, FT}, WB<:AbstractBoundary{West, FT}, + TT<:StructArray{<:TopographyElement{FT}}, } north::NB south::SB east::EB west::WB - topography::StructArray{TopographyElement{FT}} + topography::TT - function Domain{FT, NB, SB, EB, WB}( + function Domain{FT, NB, SB, EB, WB, TT}( north::NB, south::SB, east::EB, west::WB, - topography::StructArray{TopographyElement{FT}}, + topography::TT, ) where { FT<:AbstractFloat, NB<:AbstractBoundary{North, FT}, SB<:AbstractBoundary{South, FT}, EB<:AbstractBoundary{East, FT}, WB<:AbstractBoundary{West, FT}, + TT<:StructArray{<:TopographyElement{FT}}, } if !periodic_compat(north, south) throw(ArgumentError("North and south boundary walls are not \ @@ -618,7 +620,7 @@ struct Domain{ throw(ArgumentError("East boundary value is less than west \ boundary value.")) end - new{FT, NB, SB, EB, WB}(north, south, east, west, topography) + new{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) end Domain( @@ -626,15 +628,16 @@ struct Domain{ south::SB, east::EB, west::WB, - topography::StructArray{TopographyElement{FT}}, + topography::TT, ) where { FT<:AbstractFloat, NB<:AbstractBoundary{North, FT}, SB<:AbstractBoundary{South, FT}, EB<:AbstractBoundary{East, FT}, WB<:AbstractBoundary{West, FT}, + TT<:StructArray{<:TopographyElement{FT}}, } = - Domain{FT, NB, SB, EB, WB}(north, south, east, west, topography) + Domain{FT, NB, SB, EB, WB, TT}(north, south, east, west, topography) end """ @@ -659,7 +662,7 @@ Domain( EB<:AbstractBoundary{East, FT}, WB<:AbstractBoundary{West, FT}, } = - Domain{FT, NB, SB, EB, WB}( + Domain( north, south, east, @@ -869,4 +872,4 @@ Domain( (ybounds[2] - ybounds[1]) / Ny, [CellFloes{FT}() for i in 1:(Nx + 1), j in 1:(Ny + 1)], ) - \ No newline at end of file + diff --git a/src/simulation_components/model.jl b/src/simulation_components/model.jl index 0a9e8a0..6415775 100644 --- a/src/simulation_components/model.jl +++ b/src/simulation_components/model.jl @@ -45,28 +45,29 @@ simulation. This should be the length of the floes array at the beginning of the run. """ struct Model{ - FT<:AbstractFloat, - GT<:AbstractGrid{FT}, - DT<:Domain{ + FT<:AbstractFloat, # float type + GT<:AbstractGrid{FT}, # grid type + DT<:Domain{ # domain type FT, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary, <:AbstractBoundary, }, + FLT<:StructArray{<:Floe{FT}}, # floe list type } grid::GT ocean::Ocean{FT} atmos::Atmos{FT} domain::DT - floes::StructArray{Floe{FT}} # See floes.jl for floe creation + floes::FLT # See floes.jl for floe creation - function Model{FT, GT, DT}( + function Model{FT, GT, DT, FLT}( grid::GT, ocean::Ocean{FT}, atmos::Atmos{FT}, domain::DT, - floes::StructArray{Floe{FT}}, + floes::FLT, ) where { FT<:AbstractFloat, GT<:AbstractGrid{FT}, @@ -77,6 +78,7 @@ struct Model{ <:AbstractBoundary, <:AbstractBoundary, }, + FLT<:StructArray{<:Floe{FT}}, } if !domain_in_grid(domain, grid) throw(ArgumentError("Domain does not fit within grid.")) @@ -93,7 +95,7 @@ struct Model{ warmer than the ocean. This is not a situation in which the \ thermodynamics are setup for right now." end - new{FT, GT, DT}(grid, ocean, atmos, domain, floes) + new{FT, GT, DT, FLT}(grid, ocean, atmos, domain, floes) end Model( @@ -101,7 +103,7 @@ struct Model{ ocean::Ocean{FT}, atmos::Atmos{FT}, domain::DT, - floes::StructArray{Floe{FT}}, + floes::FLT, ) where { FT<:AbstractFloat, GT<:AbstractGrid{FT}, @@ -112,6 +114,7 @@ struct Model{ <:AbstractBoundary, <:AbstractBoundary, }, + FLT<:StructArray{<:Floe{FT}}, } = - Model{FT, GT, DT}(grid, ocean, atmos, domain, floes) + Model{FT, GT, DT, FLT}(grid, ocean, atmos, domain, floes) end \ No newline at end of file diff --git a/src/simulation_components/simulation.jl b/src/simulation_components/simulation.jl index c441e5b..325c82d 100644 --- a/src/simulation_components/simulation.jl +++ b/src/simulation_components/simulation.jl @@ -49,16 +49,18 @@ The user can also define settings for each physical process. """ @kwdef struct Simulation{ FT<:AbstractFloat, - GT<:AbstractGrid, - DT<:Domain, + MT<:Model{FT, <:AbstractGrid, <:Domain}, CT<:AbstractFractureCriteria, + PT<:AbstractSubFloePointsGenerator, RT<:Random.AbstractRNG, - IW<:StructVector{<:InitialStateOutputWriter}, - FW<:StructVector{<:FloeOutputWriter}, - GW<:StructVector{<:GridOutputWriter}, - CW<:StructVector{<:CheckpointOutputWriter}, + OT<:OutputWriters{ + <:StructVector{<:InitialStateOutputWriter}, + <:StructVector{<:FloeOutputWriter}, + <:StructVector{<:GridOutputWriter}, + <:StructVector{<:CheckpointOutputWriter}, + }, } - model::Model{FT, GT, DT} # Model to simulate + model::MT # Model to simulate consts::Constants{FT} = Constants() # Constants used in Simulation rng::RT = Xoshiro() # Random number generator verbose::Bool = false # String output printed during run @@ -67,12 +69,12 @@ The user can also define settings for each physical process. Δt::Int = 10 # Simulation timestep (seconds) nΔt::Int = 7500 # Total timesteps simulation runs for # Physical Processes ------------------------------------------------------- - coupling_settings::CouplingSettings = CouplingSettings() + coupling_settings::CouplingSettings{PT} = CouplingSettings() collision_settings::CollisionSettings{FT} = CollisionSettings() fracture_settings::FractureSettings{CT} = FractureSettings() simp_settings::SimplificationSettings{FT} = SimplificationSettings() # Output Writers ----------------------------------------------------------- - writers::OutputWriters{IW, FW, GW, CW} = OutputWriters() + writers::OT = OutputWriters() end """ diff --git a/test/test_floe_utils.jl b/test/test_floe_utils.jl index 7ade8fd..ce59a8f 100644 --- a/test/test_floe_utils.jl +++ b/test/test_floe_utils.jl @@ -187,12 +187,21 @@ # ------------------------- Test intersection of lines --------------------- l1 = [[[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]]] l2 = [[[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]]] - @test Subzero.intersect_lines(l1, l2) == [3.0 0.0] + @test issetequal( + Subzero.intersect_lines(l1, l2), + Set([(3.0, 0.0)]), + ) l1 = [[[0., -1], [1, 1], [2, -1], [3, 1]]] l2 = [[[0., 0], [1, 0], [3, 0]]] - @test Subzero.intersect_lines(l1, l2) == [0.5 0; 1.5 0; 2.5 0] + @test issetequal( + Subzero.intersect_lines(l1, l2), + Set([(0.5, -0.0), (1.5, 0), (2.5, -0.0)]), + ) l2 = [[[10., 10]]] - @test Subzero.intersect_lines(l1, l2) == zeros(2,0) + @test issetequal( + Subzero.intersect_lines(l1, l2), + Set{Tuple{Float64, Float64}}(), + ) # -------------- Test cutting polygon through horizontal line -------------- # Cut a hexagon through the line y = -1 diff --git a/test/test_physical_processes/test_collisions.jl b/test/test_physical_processes/test_collisions.jl index ee68cad..6f2d623 100644 --- a/test/test_physical_processes/test_collisions.jl +++ b/test/test_physical_processes/test_collisions.jl @@ -644,6 +644,7 @@ @test floe_arr[1].interactions[1, Subzero.xforce] != floe_arr[1].interactions[2, Subzero.xforce] && floe_arr[1].interactions[1, Subzero.xforce] != floe_arr[1].interactions[3,Subzero.xforce] @test floe_arr[1].interactions[1, Subzero.yforce] != floe_arr[1].interactions[2, Subzero.yforce] && floe_arr[1].interactions[1, Subzero.yforce] != floe_arr[1].interactions[3, Subzero.yforce] + empty!(small_rect.ghosts) floe_arr = StructArray([small_rect, bound_rect]) for i in eachindex(floe_arr) floe_arr.id[i] = Float64(i)