Bump actions/checkout from 6 to 7

.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 }}

.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'

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

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])
 

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()

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)