diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f3e180b..f136f68 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -26,7 +26,7 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.version }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index e5b7b7d..821c3a3 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -21,7 +21,7 @@ jobs: statuses: write pull-requests: read steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: julia-actions/setup-julia@v3 with: version: '1.12' diff --git a/src/Models/WaveGrowthModels2D.jl b/src/Models/WaveGrowthModels2D.jl index d2e9561..eb90c07 100644 --- a/src/Models/WaveGrowthModels2D.jl +++ b/src/Models/WaveGrowthModels2D.jl @@ -217,7 +217,11 @@ function WaveGrowth2D(; grid::GG, periodic_boundary=true, boundary_type="same", # or "minimal", "same", default is same, only used if periodic_boundary is false CBsets=nothing, - spline_order::Int=1, # B-spline deposition order: 1=CIC (default), 2=TSC, 3=cubic + spline_order::Int=3, # B-spline particle→grid deposition order: 1=CIC, 2=TSC, 3=cubic (default). + # Order 1/2 imprint a 4-fold grid-axis "cross" on a symmetric TC eye + # (issue #59); cubic (3) restores axisymmetry. Auto-falls back to 1 on + # tripolar grids (higher-order seam remap unsupported). Verify a pure- + # propagation ring stays circular if you change this. windsea_merge::Bool=false, # false = additive deposition (default); true = wind-sea-favoured merge! (Hanson&Phillips 2001) movie=false) where {PP<:Union{ParticleDefaults2D,String},GG<:AbstractGrid} @@ -229,9 +233,10 @@ function WaveGrowth2D(; grid::GG, error("spline_order must be 1, 2, or 3 (got $spline_order)") end # Higher-order deposition is not yet supported on tripolar grids (multi-point seam remap is a - # separate task); guard rather than silently deposit at the wrong order. + # separate task); fall back to CIC there rather than error, so the default (3) works anywhere. if spline_order > 1 && typeof(grid) <: MeshGrids && occursin("Tripolar", string(nameof(typeof(grid.stats.Ny)))) - error("spline_order > 1 is not yet supported on tripolar grids; use spline_order=1.") + @warn "spline_order > 1 is not yet supported on tripolar grids; falling back to spline_order = 1." + spline_order = 1 end # initialize state {SharedArray} given grid and layers diff --git a/src/Operators/mapping_2D.jl b/src/Operators/mapping_2D.jl index 90bc546..c9cece0 100644 --- a/src/Operators/mapping_2D.jl +++ b/src/Operators/mapping_2D.jl @@ -336,8 +336,14 @@ function remesh!(PI::ParticleInstance2D, S::StateTypeL1, # tiny, ill-defined net momentum that drives the #64 remesh limit cycle at wind/calm # interfaces, while scaling with the wind so it never over-deactivates an energetic # wind sea. α is tunable (PM is fully developed → an upper bound, so α is small). - m_amp_min = ODEs.windsea_alpha > 0 ? - max(ODEs.m_amp_minimum, ODEs.windsea_alpha * FetchRelations.windsea_momentum_PM(sqrt(wind_speed_squared))) : + # dt-scaled wind-sea floor: the momentum eroded per unit time ∝ α/dt, so a constant + # α drains the field at small dt (collapse). Scaling α ∝ min(1, dt/dt_ref) with + # dt_ref = 10 min (the validated sweet spot) bounds the erosion and prevents the + # small-dt collapse, making `windsea_alpha` safe to leave on by default. See the + # eyewall-stability FINDINGS (analysis/stability) for the α/dt collapse calibration. + α_eff = ODEs.windsea_alpha * min(1.0, ODEs.timestep / 600.0) + m_amp_min = α_eff > 0 ? + max(ODEs.m_amp_minimum, α_eff * FetchRelations.windsea_momentum_PM(sqrt(wind_speed_squared))) : ODEs.m_amp_minimum m_amp = speed(u_state[2], u_state[3]) diff --git a/src/ParticleSystems/particle_waves_v6.jl b/src/ParticleSystems/particle_waves_v6.jl index 78233d6..08c3a96 100644 --- a/src/ParticleSystems/particle_waves_v6.jl +++ b/src/ParticleSystems/particle_waves_v6.jl @@ -126,8 +126,21 @@ $(DocStringExtensions.FIELDS) wind_min_squared::Float64 = 4.0 "minimum momentum amplitude for vertex reconstruction (crash floor); below this threshold the direction is ill-defined" m_amp_minimum::Float64 = 1e-6 - "wind-sea scaling of the reconstruction floor: local threshold = max(m_amp_minimum, windsea_alpha · |m|_ws,PM(U_local)). 0 disables the local scaling (constant m_amp_minimum)." - windsea_alpha::Float64 = 0.0 + """ + Wind-sea-scaled reconstruction floor: local threshold = + max(m_amp_minimum, α_eff · |m|_ws,PM(U_local)), with α_eff = windsea_alpha · min(1, dt/10min) + automatically dt-scaled (bounds erosion at small dt; see `mapping_2D`). |m|_ws,PM is the + fully-developed Pierson–Moskowitz wind-sea momentum (∝ U³), so the floor scales with the + local wind and never over-deactivates an energetic sea; in calm regions it falls back to + `m_amp_minimum` (swell is never discarded). + + This damps the remesh limit cycle at wind/calm interfaces and the rotating TC eyewall + (energy-conserving for a resolved eyewall). DEFAULT ≈0.005 — the validated stable-TC value + (use with spline_order=3 and a resolved eyewall, dx ≲ Rmax/20, wave Courant C_cg = c_g·dt/dx + ≈ 2–5). Set 0 to disable. CAUTION: on an UNDER-resolved eyewall the floor can collapse the + storm — resolve the eyewall first. See analysis/stability FINDINGS_eyewall_stability_courant. + """ + windsea_alpha::Float64 = 0.005 "solver method for ODE system" #alternatives #Rosenbrock23(), AutoVern7(Rodas4()) ,AutoTsit5(Rosenbrock23()) , Tsit5() diff --git a/test/manual/T04_2D_growing_decaying_winds.jl b/test/manual/T04_2D_growing_decaying_winds.jl index c6a9b79..6fc8ce4 100644 --- a/test/manual/T04_2D_growing_decaying_winds.jl +++ b/test/manual/T04_2D_growing_decaying_winds.jl @@ -37,7 +37,7 @@ using PiCLES.Plotting.movie: init_movie_2D_box_plot #sign.(rand(-1:1, 10, 10)) -save_path = "plots/tests/T04_2D_growing_decaying_winds/" +save_path = "PiCLES/plots/tests/T04_2D_growing_decaying_winds_P2/" mkpath(save_path) ##### basic parameters @@ -92,7 +92,7 @@ ODE_settings = PW.ODESettings( function make_reg_test(wave_model, save_path; plot_name="dummy", N=36, axline=0) ### build Simulation - wave_simulation = Simulation(wave_model, Δt=DT, stop_time=1hour)#1hours) + wave_simulation = Simulation(wave_model, Δt=DT, stop_time=5hour)#1hours) initialize_simulation!(wave_simulation) # run simulation @@ -116,7 +116,7 @@ end # % half domain tests -gridmesh = [(i, j) for i in [-10,10], j in [0]] +gridmesh = [(i, j) for i in [-30,-20,-10,10], j in [0]] #gridmesh = [(i, j) for i in [10], j in [0]] #for I in CartesianIndices(gridmesh) @@ -149,9 +149,10 @@ for (U10, V10) in gridmesh boundary_type="same", minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), # # minimal_state=FetchRelations.MinimalState(2, 2, DT) * 1, + spline_order=2, movie=true) - make_reg_test(wave_model, save_path, plot_name="T02_2D_growing_U" * string(U10) * "_V" * string(V10), N=40, axline=x0/1e3) + make_reg_test(wave_model, save_path, plot_name="T02_2D_growing_U" * string(U10) * "_V" * string(V10), N=48*6, axline=x0/1e3) end # %% @@ -186,6 +187,7 @@ for (U10, V10) in gridmesh boundary_type="same", minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), # # minimal_state=FetchRelations.MinimalState(2, 2, DT) * 1, + spline_order=2, movie=true) make_reg_test(wave_model, save_path, plot_name="T02_2D_decaying_U" * string(U10) * "_V" * string(V10), N=60, axline=x0/1e3) @@ -223,6 +225,7 @@ for (U10, V10) in gridmesh boundary_type="same", minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), # # minimal_state=FetchRelations.MinimalState(2, 2, DT) * 1, + spline_order=2, movie=true) make_reg_test(wave_model, save_path, plot_name="T02_2D_divergence_U" * string(U10) * "_V" * string(V10), N=100, axline=x0/1e3)