diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 6f52ed5..f82af94 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -1,7 +1,7 @@ name: CompatHelper on: schedule: - - cron: 0 0 * * * + - cron: 0 14 * * * workflow_dispatch: permissions: contents: write diff --git a/Project.toml b/Project.toml index 17c5c26..662b9f8 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ 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,6 +13,7 @@ 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" diff --git a/examples/test_run.jl b/examples/test_run.jl index 6cfe8f8..a40cca7 100644 --- a/examples/test_run.jl +++ b/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() diff --git a/src/physical_processes/coupling.jl b/src/physical_processes/coupling.jl index d311c30..3e48a62 100644 --- a/src/physical_processes/coupling.jl +++ b/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 --------------# diff --git a/src/tools/plotting.jl b/src/tools/plotting.jl index 21ff51c..ae94982 100644 --- a/src/tools/plotting.jl +++ b/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]), diff --git a/test/test_floe.jl b/test/test_floe.jl index 3dd22b9..11c5dd7 100644 --- a/test/test_floe.jl +++ b/test/test_floe.jl @@ -6,51 +6,7 @@ close(file) poly1 = LG.Polygon(Subzero.valid_polyvec!(floe_coords[1])) centroid1 = LG.GeoInterface.coordinates(LG.centroid(poly1)) - origin_coords = Subzero.translate( - floe_coords[1], - -centroid1[1], - -centroid1[2], - ) - xo, yo = Subzero.separate_xy(origin_coords) - rmax = sqrt(maximum([sum(xo[i]^2 + yo[i]^2) for i in eachindex(xo)])) area = LG.area(poly1) - mc_x, mc_y, status = Subzero.generate_subfloe_points( - MonteCarloPointsGenerator(), - origin_coords, - rmax, - area, - Subzero.Status(), - Xoshiro(1) - ) - @test length(mc_x) == length(mc_y) && length(mc_x) > 0 - in_on = inpoly2(hcat(mc_x, mc_y), hcat(xo, yo)) - mc_in = in_on[:, 1] .| in_on[:, 2] - @test all(mc_in) - @test abs(sum(mc_in)/1000 * 4 * rmax^2 - area)/area < 0.1 - @test status.tag == Subzero.active - # Test that random number generator is working - mc_x2, mc_y2, status2 = Subzero.generate_subfloe_points( - MonteCarloPointsGenerator(), - origin_coords, - rmax, - area, - Subzero.Status(), - Xoshiro(1) - ) - @test all(mc_x .== mc_x2) - @test all(mc_y .== mc_y2) - @test status2.tag == Subzero.active - - mc_x3, mc_y3, status3 = Subzero.generate_subfloe_points( - MonteCarloPointsGenerator{Float32}(), - origin_coords, - rmax, - area, - Subzero.Status(), - Xoshiro(1) - ) - @test status3.tag == Subzero.active - @test eltype(mc_x3) == eltype(mc_y3) == Float32 # Test InteractionFields enum interactions = range(1, 7)' diff --git a/test/test_physical_processes/test_coupling.jl b/test/test_physical_processes/test_coupling.jl index 6fc8a0a..e710165 100644 --- a/test/test_physical_processes/test_coupling.jl +++ b/test/test_physical_processes/test_coupling.jl @@ -1,5 +1,160 @@ @testset "Coupling" begin FT = Float64 + @testset "Generating sub-floe points" begin + FT = Float64 + # test generating monte carlo points + file = jldopen("inputs/floe_shapes.jld2", "r") + floe_coords = file["floe_vertices"][1:end] + close(file) + poly1 = LG.Polygon(Subzero.valid_polyvec!(floe_coords[1])) + centroid1 = LG.GeoInterface.coordinates(LG.centroid(poly1)) + origin_coords = Subzero.translate( + floe_coords[1], + -centroid1[1], + -centroid1[2], + ) + xo, yo = Subzero.separate_xy(origin_coords) + rmax = sqrt(maximum([sum(xo[i]^2 + yo[i]^2) for i in eachindex(xo)])) + area = LG.area(poly1) + mc_x, mc_y, status = Subzero.generate_subfloe_points( + MonteCarloPointsGenerator(), + origin_coords, + rmax, + area, + Subzero.Status(), + Xoshiro(1) + ) + @test length(mc_x) == length(mc_y) && length(mc_x) > 0 + in_on = inpoly2(hcat(mc_x, mc_y), hcat(xo, yo)) + mc_in = in_on[:, 1] .| in_on[:, 2] + @test all(mc_in) + @test abs(sum(mc_in)/1000 * 4 * rmax^2 - area)/area < 0.1 + @test status.tag == Subzero.active + # Test that random number generator is working + mc_x2, mc_y2, status2 = Subzero.generate_subfloe_points( + MonteCarloPointsGenerator(), + origin_coords, + rmax, + area, + Subzero.Status(), + Xoshiro(1) + ) + @test all(mc_x .== mc_x2) + @test all(mc_y .== mc_y2) + @test status2.tag == Subzero.active + + mc_x3, mc_y3, status3 = Subzero.generate_subfloe_points( + MonteCarloPointsGenerator{Float32}(), + origin_coords, + rmax, + area, + Subzero.Status(), + Xoshiro(1) + ) + @test status3.tag == Subzero.active + @test eltype(mc_x3) == eltype(mc_y3) == Float32 + + # test generating sub-grid points for grid with Δx = Δy = 10 + point_generator = SubGridPointsGenerator{Float64}(10/sqrt(2)) + # Floe is smaller than grid cells --> centroid and vertices added + square = [[ + [-2.5, -2.5], + [-2.5, 2.5], + [2.5, 2.5], + [2.5, -2.5], + [-2.5, -2.5], + ]] + xpoints, ypoints = Subzero.generate_subfloe_points( + point_generator, + square, + 0.0, # Not used + 0.0, # Not used + Subzero.Status(), + Xoshiro(), # Not random + ) + @test xpoints == [-2.5, -2.5, 2.5, 2.5, 0.0] + @test ypoints == [-2.5, 2.5, 2.5, -2.5, 0.0] + # Floe is larger than grid cell + tall_rect = [[ + [-2.0, -10.0], + [-2.0, 10.0], + [2.0, 10.0], + [2.0, -10.0], + [-2.0, -10.0], + ]] + xpoints, ypoints = Subzero.generate_subfloe_points( + point_generator, + tall_rect, + 0.0, # Not used + 0.0, # Not used + Subzero.Status(), + Xoshiro(), # Not random + ) + @test xpoints == [repeat([-2.0], 5); repeat([2.0], 5); repeat([0.0], 3)] + @test all(isapprox.( + ypoints, + [-10.0, -6.46447, 0.0, 6.46447, 10.0, 10.0, 6.46447, 0.0, + -6.46447, -10, -6.46447, 0.0, 6.46447, + ], + atol = 1e-5)) + + wide_rect = [[ + [-10.0, -2.0], + [-10.0, 2.0], + [10.0, 2.0], + [10.0, -2.0], + [-10.0, -2.0], + ]] + xpoints, ypoints = Subzero.generate_subfloe_points( + point_generator, + wide_rect, + 0.0, # Not used + 0.0, # Not used + Subzero.Status(), + Xoshiro(), # Not random + ) + @test all(isapprox.( + xpoints, + [-10, -10, -6.46447, 0.0, 6.46447, 10, 10, 6.46447, 0.0, + -6.464466, -6.46447, 0, 6.46447, + ], + atol = 1e-5 + )) + @test ypoints == [-2; repeat([2], 5); repeat([-2], 4); repeat([0], 3)] + trapeziod = [[ + [-8.0, -8.0], + [-4.0, 8.0], + [4.0, 8.0], + [8.0, -8.0], + [-8.0, -8.0], + ]] + xpoints, ypoints = Subzero.generate_subfloe_points( + point_generator, + trapeziod, + 0.0, # Not used + 0.0, # Not used + Subzero.Status(), + Xoshiro(), # Not random + ) + @test all(isapprox.( + xpoints, + [-8, -7.14251, -6.0, -4.85749, -4.0, 0.0, 4.0, 4.85749, 6.0, + 7.14251, 8.0, 4.46447, 0.0, -4.46447, -4.46447, 0.0, 4.46447, + -4.46447, 0.0, 4.46447, -4.46447, 0.0, 4.46447 + ], + atol = 1e-5 + )) + @test all(isapprox.( + ypoints, + [-8; -4.57003; 0.0; 4.57003; repeat([8.0], 3); 4.57003; 0.0; + -4.57003; repeat([-8.0], 4); repeat([-4.46447], 3); + repeat([0.0], 3); repeat([4.46447], 3); + ], + atol = 1e-5 + )) + + end + @testset "Coupling Helper Functions" begin grid = Subzero.RegRectilinearGrid( (-10, 10),