Merge pull request #3379 from CliMA/gb/restart2

Warn when dt_rad and dt_checkpoints are not compatible

src/callbacks/get_callbacks.jl

@@ -209,14 +209,17 @@ function get_callbacks(config, sim_info, atmos, params, Y, p, t_start)
         ),
     )
 
+    # Save dt_save_state_to_disk as a Dates.Period object. This is used to check
+    # if it is an integer multiple of other frequencies.
+    dt_save_state_to_disk_dates = Dates.today() # Value will be overwritten
     if occursin("months", parsed_args["dt_save_state_to_disk"])
         months = match(r"^(\d+)months$", parsed_args["dt_save_state_to_disk"])
         isnothing(months) && error(
             "$(period_str) has to be of the form <NUM>months, e.g. 2months for 2 months",
         )
-        period_dates = Dates.Month(parse(Int, first(months)))
+        dt_save_state_to_disk_dates = Dates.Month(parse(Int, first(months)))
         schedule = CAD.EveryCalendarDtSchedule(
-            period_dates;
+            dt_save_state_to_disk_dates;
             reference_date = p.start_date,
             date_last = p.start_date + Dates.Second(t_start),
         )
@@ -231,6 +234,9 @@ function get_callbacks(config, sim_info, atmos, params, Y, p, t_start)
         dt_save_state_to_disk =
             time_to_seconds(parsed_args["dt_save_state_to_disk"])
         if !(dt_save_state_to_disk == Inf)
+            # We use Millisecond to support fractional seconds, eg. 0.1
+            dt_save_state_to_disk_dates =
+                Dates.Millisecond(dt_save_state_to_disk)
             callbacks = (
                 callbacks...,
                 call_every_dt(
@@ -273,6 +279,14 @@ function get_callbacks(config, sim_info, atmos, params, Y, p, t_start)
 
     if atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode
         dt_rad = FT(time_to_seconds(parsed_args["dt_rad"]))
+        # We use Millisecond to support fractional seconds, eg. 0.1
+        dt_rad_ms = Dates.Millisecond(dt_rad)
+        if parsed_args["dt_save_state_to_disk"] != "Inf" &&
+           !CA.isdivisible(dt_save_state_to_disk_dates, dt_rad_ms)
+            @warn "Radiation period ($(dt_rad_ms)) is not an even divisor of the checkpoint frequency ($dt_save_state_to_disk_dates)"
+            @warn "This simulation will not be reproducible when restarted"
+        end
+
         callbacks =
             (callbacks..., call_every_dt(rrtmgp_model_callback!, dt_rad))
     end

src/utils/utilities.jl

@@ -403,3 +403,54 @@ function gaussian_smooth(arr::AbstractArray, sigma::Int = 1)
 
     return smoothed_arr
 end
+
+"""
+    isdivisible(dt_large::Dates.Period, dt_small::Dates.Period)
+
+Check if two periods are evenly divisible, i.e., if the larger period can be
+expressed as an integer multiple of the smaller period.
+
+In this, take into account the case when periods do not have fixed size, e.g.,
+one month is a variable number of days.
+
+# Examples
+```
+julia> isdivisible(Dates.Year(1), Dates.Month(1))
+true
+
+julia> isdivisible(Dates.Month(1), Dates.Day(1))
+true
+
+julia> isdivisible(Dates.Month(1), Dates.Week(1))
+false
+```
+
+## Notes
+
+Not all the combinations are fully implemented. If something is missing, please
+consider adding it.
+"""
+function isdivisible(dt_large::Dates.Period, dt_small::Dates.Period)
+    @warn "The combination $(typeof(dt_large)) and $(dt_small) was not covered. Please add a method to handle this case."
+    return false
+end
+
+# For FixedPeriod and OtherPeriod, it is easy, we can directly divide the two
+# (as long as they are both the same)
+function isdivisible(dt_large::Dates.FixedPeriod, dt_small::Dates.FixedPeriod)
+    return isinteger(dt_large / dt_small)
+end
+
+function isdivisible(dt_large::Dates.OtherPeriod, dt_small::Dates.OtherPeriod)
+    return isinteger(dt_large / dt_small)
+end
+
+function isdivisible(
+    dt_large::Union{Dates.Month, Dates.Year},
+    dt_small::Dates.FixedPeriod,
+)
+    # The only case where periods are commensurate for Month/Year is when we
+    # have a Day or an integer divisor of a day. (Note that 365 and 366 don't
+    # have any common divisor)
+    return isinteger(Dates.Day(1) / dt_small)
+end

test/restart.jl

@@ -241,7 +241,7 @@ function test_restart(test_dict; job_id, comms_ctx, more_ignore = Symbol[])
 
     simulation_restarted = CA.get_simulation(config_should_be_same)
 
-    if pkgversion(RRTMGP) <= v"0.19"
+    if pkgversion(RRTMGP) < v"0.20"
         # Versions of RRTMGP older than 0.20 have a bug and do not set the
         # flux_dn_dir, so that face_clear_sw_direct_flux_dn and
         # face_sw_direct_flux_dn is uninitialized and not deterministic

test/utilities.jl

@@ -1,6 +1,7 @@
 using Test
 using ClimaComms
 ClimaComms.@import_required_backends
+import Dates
 using Random
 Random.seed!(1234)
 import ClimaAtmos as CA
@@ -30,6 +31,18 @@ end
     @test extrema(randy)[2] >= smoothed[2]
 end
 
+@testset "isdivisible" begin
+    @test CA.isdivisible(Dates.Month(1), Dates.Day(1))
+    @test !CA.isdivisible(Dates.Month(1), Dates.Day(25))
+    @test CA.isdivisible(Dates.Week(1), Dates.Day(1))
+    @test CA.isdivisible(Dates.Day(1), Dates.Hour(1))
+    @test CA.isdivisible(Dates.Hour(1), Dates.Second(1))
+    @test CA.isdivisible(Dates.Minute(1), Dates.Second(30))
+    @test !CA.isdivisible(Dates.Minute(1), Dates.Second(13))
+    @test !CA.isdivisible(Dates.Day(1), Dates.Second(1e6))
+    @test CA.isdivisible(Dates.Month(1), Dates.Hour(1))
+end
+
 @testset "kinetic_energy (c.f. analytical function)" begin
     # Test kinetic energy function for staggered grids
     # given an analytical expression for the velocity profiles