Debug subfloe points (#77)

* Save ocean tau

* Fix sub-grid algorithm

* Update subgrid algorithm

* Basic test of new subgrid point algorithm

* Finish new subgrid algorithm

* Update bad test

.github/workflows/CompatHelper.yml

@@ -1,7 +1,7 @@
 name: CompatHelper
 on:
  schedule:
- - cron: 0 0 * * *
+ - cron: 0 14 * * *
  workflow_dispatch:
 permissions:
  contents: write

Project.toml

@@ -6,14 +6,14 @@ 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"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
+LibGEOSMakie = "1902729b-b4d4-4bb0-ae7d-e2f19e2a15fb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -23,7 +23,6 @@ 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"

examples/test_run.jl

@@ -1,54 +1,11 @@
 using JLD2, Random, Statistics, Subzero, BenchmarkTools, StructArrays
 import LibGEOS as LG
 
-Δt = 10
-max_overlap = 0.75
-Lx = 1e5
-Ly = Lx
-hmean = 0.25
-Δh = 0.0
-grid = RegRectilinearGrid(
- (-Lx, Lx),
- (-Ly, Ly),
- 1e4,
- 1e4,
-)
-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]),
-)
-# Diagonal floe barely overlaping with eastern collision boundary
-efloe_small = Floe(
- [[
- [9.5e4, 0.0],
- [9e4, 0.5e4],
- [10e4, 2.5e4],
- [10.05e4, 2e4],
- [9.5e4, 0.0],
- ]],
- hmean,
- Δh,
-)
-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
 const Lx = 1e5
 const Ly = 1e5
-const Δgrid = 1e4
+const Δgrid = 2e3
 const hmean = 0.25
 const Δh = 0.0
 const Δt = 10
@@ -59,74 +16,103 @@ grid = RegRectilinearGrid(
  Δgrid,
  Δgrid,
 )
-zero_ocn = Ocean(grid, 0.0, 0.0, 0.0)
+zonal_ocn = Ocean(grid, 0.5, 0.0, 0.0)
 
-zero_atmos = Atmos(grid, 0.0, 0.0, 0.0)
+zero_atmos = Atmos(grid, 0.5, 0.0, 0.0)
 
 
 domain = Subzero.Domain(
- CollisionBoundary(North, grid),
- CollisionBoundary(South, grid),
- CompressionBoundary(East, grid, -0.1),
- CompressionBoundary(West, grid, 0.1),
+ PeriodicBoundary(North, grid),
+ PeriodicBoundary(South, grid),
+ PeriodicBoundary(East, grid),
+ PeriodicBoundary(West, grid),
 )
 
 # Floe instantiation
-coords = [[[0.0, 0.0], [0.0, 1e5], [1e5, 1e5], [1e5, 0.0], [0.0, 0.0]]]
+coords = [[
+ [0.0, 0.0],
+ [5e4, 1e5],
+ [1e5, 1e5],
+ [1e5, 0.0],
+ [0.0, 0.0],
+]]
+
+coupling_settings = CouplingSettings(
+ subfloe_point_generator = SubGridPointsGenerator(grid, 1),
+ two_way_coupling_on = true,
+)
 
 floe_arr = initialize_floe_field(
  FT,
- 50,
- [1.0],
+ 30,
+ [0.8],
  domain,
  0.5,
- 0.0,
- min_floe_area = 1e6,
+ 0.0;
+ coupling_settings = coupling_settings,
+ rng = Xoshiro(1),
 )
-
+floe_arr.α .= 1e-5
 model = Model(
  grid,
- zero_ocn,
+ zonal_ocn,
  zero_atmos,
  domain,
  floe_arr,
 )
 dir = "output/sim"
 writers = OutputWriters(
- InitialStateOutputWriter(
- dir = dir,
- overwrite = true
- ),
- FloeOutputWriter(
- 250,
- dir = dir,
- overwrite = true,
- ),
+ # InitialStateOutputWriter(
+ #  dir = dir,
+ #  overwrite = true
+ # ),
+ # FloeOutputWriter(
+ #  250,
+ #  dir = dir,
+ #  overwrite = true,
+ # ),
 )
 simulation = Simulation(
  name = "sim",
  model = model,
  Δt = 10,
- nΔt = 4000,
+ nΔt = 200,
  writers = writers,
  verbose = true,
- consts = Constants(Cd_io = 0, Cd_ia = 0, Cd_ao = 0),
- fracture_settings = FractureSettings(
- fractures_on = true,
- criteria = HiblerYieldCurve(model.floes),
- Δt = 75,
- )
+ consts = Constants(Cd_ia = 0, Cd_ao = 0),
+ coupling_settings = coupling_settings,
+
 )
 
 #@benchmark timestep_sim!(simulation, 10) setup=(sim=deepcopy(simulation))
+function run_with_saving(simulation, tstep_between_save)
+ nsave = fld(simulation.nΔt, tstep_between_save) + 1
+ τx = zeros(FT, grid.Nx + 1, grid.Ny + 1, nsave)
+ Subzero.startup_sim(simulation, nothing, 1)
+ tstep = 0
+ counter = 1
+ while tstep <= simulation.nΔt
+ # Timestep the simulation forward
+ timestep_sim!(simulation, tstep)
+ if mod(tstep, tstep_between_save) == 0
+ τx[:, :, counter] .= simulation.model.ocean.τx
+ counter += 1
+ end
+ tstep+=1
+ end
+ Subzero.teardown_sim(simulation)
+ return τx
+end
+
+τx = run_with_saving(simulation, 50)
 
 
-time_run(simulation) = @time run!(simulation)
-time_run(simulation)
+# time_run(simulation) = @time run!(simulation)
+# time_run(simulation)
 
-Subzero.create_sim_gif("output/sim/floes.jld2", 
- "output/sim/initial_state.jld2",
- "output/sim/sim.gif")
+# Subzero.create_sim_gif("output/sim/floes.jld2", 
+#  "output/sim/initial_state.jld2",
+#  "output/sim/sim.gif")
 # # Run simulation
 #time_run(simulation)
 #Profile.Allocs.clear()

src/physical_processes/coupling.jl

@@ -139,7 +139,9 @@ SubGridPointsGenerator{FT}(
  grid::RegRectilinearGrid,
  npoint_per_cell::Int,
 ) where {FT <: AbstractFloat} =
- SubGridPointsGenerator{FT}(min(grid.Δx, grid.Δy) / npoint_per_cell)
+ SubGridPointsGenerator{FT}(
+ min(grid.Δx, grid.Δy) / npoint_per_cell / sqrt(2)
+ )
 
 """
  generate_subfloe_points(
@@ -234,51 +236,95 @@ function generate_subfloe_points(
  status,
  rng
 ) where {FT <: AbstractFloat}
- xmax = FT(-Inf)
- ymax = FT(-Inf)
- xmin = FT(Inf)
- ymin = FT(Inf)
- @views for vert in coords[1]
- if vert[1] > xmax
- xmax = vert[1]
- end
- if vert[1] < xmin
- xmin = vert[1]
+ xmax = coords[1][1][1]
+ xmin = coords[1][1][1]
+ ymax = coords[1][1][2]
+ ymin = coords[1][1][2]
+ nverts = length(coords[1])
+
+ xpoints = Vector{FT}()
+ ypoints = Vector{FT}()
+ # Add points along edges
+ for i in 1:(nverts - 1)
+ # Determine points on edges
+ x1, y1 = coords[1][i]
+ x2, y2 = coords[1][i+1]
+ Δx = x2 - x1
+ Δy = y2 - y1
+ l = sqrt((Δx)^2 + (Δy)^2)
+ # Find maximum and minimum x and y points
+ if x1 > xmax
+ xmax = x1
+ elseif x1 < xmin
+ xmin = x1
  end
- if vert[2] > ymax
- ymax = vert[2]
+ if y1 > ymax
+ ymax = y1
+ elseif y1 < ymin
+ ymin = y1
  end
- if vert[2] < ymin
- ymin = vert[2]
+ # Add current vertex
+ push!(xpoints, x1)
+ push!(ypoints, y1)
+ # If distance between i and i+1 vertex is less than 2 sub-grid cells
+ # but greater than one, add a point inbetween those two vertices
+ if l <= 2point_generator.Δg
+ if l > point_generator.Δg
+ push!(xpoints, x1 + Δx/2)
+ push!(ypoints, y1 + Δy/2)
+ end
+ else # The edge needs more points than the corners and midpoint
+ if Δx == 0
+ y1 += point_generator.Δg/2 * sign(Δy)
+ y2 -= point_generator.Δg/2 * sign(Δy)
+ elseif Δy == 0
+ x1 += point_generator.Δg/2 * sign(Δx)
+ x2 -= point_generator.Δg/2 * sign(Δx)
+ else # shift points to still be on the line
+ m = Δy / Δx
+ x_shift = sqrt(point_generator.Δg^2 / 4(1 + m^2))
+ y_shift = m * x_shift
+ x1 += x_shift
+ x2 -= x_shift
+ y1 += y_shift
+ y2 -= y_shift
+ end
+ l = sqrt((x2 - x1)^2 + (y2 - y1)^2)
+ n_edge_points = ceil(Int, l / point_generator.Δg) + 1
+ append!(xpoints, range(x1, x2, length = n_edge_points))
+ append!(ypoints, range(y1, y2, length = n_edge_points))
  end
  end
-
- n_xpoints = round(Int, (xmax - xmin) / point_generator.Δg) + 1
- n_ypoints = round(Int, (ymax - ymin) / point_generator.Δg) + 1
- xpoints = if n_xpoints < 3
+ # Add points in the interior of the floe
+ n_xpoints = ceil(Int, (xmax - xmin) / point_generator.Δg)
+ n_ypoints = ceil(Int, (ymax - ymin) / point_generator.Δg)
+ x_interior_points = if n_xpoints < 3
+ n_xpoints = 1
  FT(0):FT(0) # coords are centered at the origin
  else
- range(xmin, xmax, length = n_xpoints)
+ range(
+ xmin + point_generator.Δg/2,
+ xmax - point_generator.Δg/2,
+ length = n_xpoints,
+ )
  end
- ypoints = if n_ypoints < 3
- FT(0):FT(0) # coords are centered at the origin
+ y_interior_points = if n_ypoints < 3
+ n_ypoints = 1
+ FT(0):FT(0)
  else
- range(ymin, ymax, length = n_ypoints)
+ range(
+ ymin + point_generator.Δg/2,
+ ymax - point_generator.Δg/2,
+ length = n_ypoints,
+ )
  end
- x_sub_floe = repeat(xpoints, length(ypoints))
- y_sub_floe = repeat(ypoints, inner = length(xpoints))
+ x_sub_floe = repeat(x_interior_points, n_ypoints)
+ y_sub_floe = repeat(y_interior_points, inner = n_xpoints)
  in_floe = points_in_poly(hcat(x_sub_floe, y_sub_floe), coords)
 
- x_sub_floe = x_sub_floe[in_floe]
- y_sub_floe = y_sub_floe[in_floe]
- if length(x_sub_floe) < 2
- x_sub_floe = FT(0):FT(0) # coords are centered at the origin
- end
- if length(y_sub_floe) < 2
- y_sub_floe = FT(0):FT(0) # coords are centered at the origin
- end
-
- return x_sub_floe, y_sub_floe, status
+ append!(xpoints, x_sub_floe[in_floe])
+ append!(ypoints, y_sub_floe[in_floe])
+ return xpoints, ypoints, status
 end
 
 #-------------- Monte Carlo Point Calculations --------------#

src/tools/plotting.jl

@@ -268,7 +268,7 @@ function create_coupled_ro_sim_gif(
  plt,
  xc/1000,
  yc/1000,
- Ro, # should be transposed 
+ Ro', # should be transposed 
  xlabel="x [km]",
  ylabel="y [km]",
  title=prettytime(t[i]),