Merge #105

105: Daily interpolation from monthly data (v2) r=LenkaNovak a=LenkaNovak

# PULL REQUEST
- Update from #91

## Purpose and Content
Adds capability to interpolate linearly between two monthly `Fields`. The demo uses prerequisite PRs outlined in #88 to obtain monthly data. 

## Benefits and Risks
- benefit: last piece for enabling specifying temporarily varying boundary conditions from a file.
- risk: performance slow-down - check are in place (see below)

## Linked Issues
- SDI: #74 
- Design: #88

## Dependent PRs
- no not merge before  #87 ,  #89 and #90

## Performance of the final solution 
- time loop only over added infrastructure of #88 : adding 0.24s per simulation day at `dt = 200s`

## PR Checklist
- [x] This PR has a corresponding issue OR is linked to an SDI.
- [x] I have followed CliMA's codebase [contribution](https://clima.github.io/ClimateMachine.jl/latest/Contributing/) and [style](https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/) guidelines OR N/A.
- [x] I have followed CliMA's [documentation policy](https://github.com/CliMA/policies/wiki/Documentation-Policy).
- [x] I have checked all issues and PRs and I certify that this PR does not duplicate an open PR.
- [x] I linted my code on my local machine prior to submission OR N/A.
- [x] Unit tests are included OR N/A.
- [x] Code used in an integration test OR N/A.
- [x] All tests ran successfully on my local machine OR N/A.
- [x] All classes, modules, and function contain docstrings OR N/A.
- [x] Documentation has been added/updated OR N/A.

Co-authored-by: LenkaNovak <[email protected]>

experiments/AMIP/moist_mpi_earth/coupler_driver.jl

@@ -87,13 +87,13 @@ if mode_name == "amip"
         sst_data,
         "SST",
         boundary_space,
-        interpolate_daily = false,
+        interpolate_daily = true,
         scaling_function = clean_sst,
         land_mask = land_mask,
         date0 = date0,
     )
     update_midmonth_data!(date0, SST_info)
-    SST_init = deepcopy(SST_info.monthly_fields[1])
+    SST_init = interpolate_midmonth_to_daily(date0, SST_info)
     ocean_params = OceanSlabParameters(FT(20), FT(1500.0), FT(800.0), FT(280.0), FT(1e-3), FT(1e-5), FT(0.06))
     ocean_sim = (; integrator = (; u = (; T_sfc = SST_init), p = (; params = ocean_params), SST_info = SST_info))
 
@@ -109,7 +109,7 @@ if mode_name == "amip"
         date0 = date0,
     )
     update_midmonth_data!(date0, SIC_info)
-    SIC_init = deepcopy(SIC_info.monthly_fields[1])
+    SIC_init = interpolate_midmonth_to_daily(date0, SIC_info)
     ice_mask = get_ice_mask.(SIC_init, FT) # here 50% and lower is considered ice free
     ice_sim = ice_init(FT; tspan = tspan, dt = Δt_cpl, space = boundary_space, saveat = saveat, ice_mask = ice_mask)
     mode_specifics = (; name = mode_name, SST_info = SST_info, SIC_info = SIC_info)
@@ -186,11 +186,11 @@ function solve_coupler!(cs, energy_check)
             @calendar_callback :(update_midmonth_data!(cs.dates.date[1], cs.mode.SST_info)) cs.dates.date[1] next_date_in_file(
                 cs.mode.SST_info,
             )
-            SST = ocean_sim.integrator.u.T_sfc .= cs.mode.SST_info.monthly_fields[1]
+            SST = ocean_sim.integrator.u.T_sfc .= interpolate_midmonth_to_daily(cs.dates.date[1], cs.mode.SST_info)
             @calendar_callback :(update_midmonth_data!(cs.dates.date[1], cs.mode.SIC_info)) cs.dates.date[1] next_date_in_file(
                 cs.mode.SIC_info,
             )
-            SIC = cs.mode.SIC_info.monthly_fields[1]
+            SIC = interpolate_midmonth_to_daily(cs.dates.date[1], cs.mode.SIC_info)
             ice_mask = ice_sim.integrator.p.Ya.ice_mask .= get_ice_mask.(SIC, FT)
 
         end
@@ -263,4 +263,3 @@ rm(REGRID_DIR; recursive = true, force = true)
 # - cs needs to be global for the monthly macro - explote other solutions
 # - SST_init is modified with SST_info even with deepcopy...
 # - replace if statements with dipatches, write better abstractions
-# - unit test for monthly file update 

experiments/AMIP/moist_mpi_earth/coupler_utils/bcfile_reader.jl

@@ -17,10 +17,10 @@ The inputs are:
     monthly_fields::C               # Tuple of the two monthly fields, that will be used for the daily interpolation
     segment_length::Vector{Int}     # length of each month segment (used in the daily interpolation)
     scaling_function::O             # function that scales, offsets or transforms the raw variable
-    interpolate_daily::I            # switch to trigger daily interpolation
-    land_mask::M                     # mask with 1 = land, 0 = ocean / sea-ice
+    interpolate_daily::Bool         # switch to trigger daily interpolation
+    land_mask::M                    # mask with 1 = land, 0 = ocean / sea-ice
 """
-struct BCFileInfo{F, S, V, D, C, O, I, M}
+struct BCFileInfo{F, S, V, D, C, O, M}
     FT::F
     datafile_cgll::S
     varname::V
@@ -31,7 +31,7 @@ struct BCFileInfo{F, S, V, D, C, O, I, M}
     monthly_fields::C
     segment_length::Vector{Int}
     scaling_function::O
-    interpolate_daily::I
+    interpolate_daily::Bool
     land_mask::M
 end
 
@@ -83,12 +83,12 @@ function bcfile_info_init(
         varname,
         weightfile,
         data_dates,
-        segment_idx0 .- Int(1),
+        deepcopy(segment_idx0),
         segment_idx0,
         current_fields,
         segment_length,
         scaling_function,
-        [interpolate_daily],
+        interpolate_daily,
         land_mask,
     )
 
@@ -105,7 +105,6 @@ Extracts boundary condition data from regridded (to model grid) NetCDF files (wh
 function update_midmonth_data!(date, bcf_info)
 
     # monthly count
-    bcf_info.segment_idx[1] += Int(1)
 
     all_dates = bcf_info.all_dates
     midmonth_idx = bcf_info.segment_idx[1]
@@ -133,7 +132,10 @@ function update_midmonth_data!(date, bcf_info)
                 ),
             Tuple(1:length(monthly_fields)),
         )
+        bcf_info.segment_length .= FT(0)
     elseif (midmonth_idx == midmonth_idx0) && (Dates.days(date - all_dates[midmonth_idx]) < 0) # for init
+        midmonth_idx = bcf_info.segment_idx[1] -= Int(1)
+        midmonth_idx = midmonth_idx < Int(1) ? midmonth_idx + Int(1) : midmonth_idx
         @warn "this time period is before BC data - using file from $(all_dates[midmonth_idx0])"
         map(
             x ->
@@ -149,6 +151,7 @@ function update_midmonth_data!(date, bcf_info)
                 ),
             Tuple(1:length(monthly_fields)),
         )
+        bcf_info.segment_length .= FT(0)
     elseif Dates.days(date - all_dates[end - 1]) > 0 # for fini
         @warn "this time period is after BC data - using file from $(all_dates[end - 1])"
         map(
@@ -165,13 +168,15 @@ function update_midmonth_data!(date, bcf_info)
                 ),
             Tuple(1:length(monthly_fields)),
         )
-    elseif Dates.days(date - all_dates[Int(midmonth_idx)]) > 2 # throw error when there are closer initial indices for the bc file data that matches this date0
+        bcf_info.segment_length .= FT(0)
+    elseif Dates.days(date - all_dates[Int(midmonth_idx + 1)]) > 2 # throw error when there are closer initial indices for the bc file data that matches this date0
         nearest_idx =
             argmin(abs.(parse(FT, datetime_to_strdate(date)) .- parse.(FT, datetime_to_strdate.(all_dates[:]))))
         @error "init data does not correspond to start date. Try initializing with `SIC_info.segment_idx = midmonth_idx = midmonth_idx0 = $nearest_idx` for this start date" # TODO: do this automatically w a warning
-    elseif Dates.days(date - all_dates[Int(midmonth_idx - 1)]) > 0 # date crosses to the next month
-        @warn "updating monthly data file"
-        bcf_info.segment_length .= Dates.days(all_dates[Int(midmonth_idx + 1)] - all_dates[Int(midmonth_idx)])
+    elseif Dates.days(date - all_dates[Int(midmonth_idx)]) > 0 # date crosses to the next month
+        midmonth_idx = bcf_info.segment_idx[1] += Int(1)
+        @warn "On $date updating monthly data files: mid-month dates = [ $(all_dates[Int(midmonth_idx)]) , $(all_dates[Int(midmonth_idx+1)]) ]"
+        bcf_info.segment_length .= (all_dates[Int(midmonth_idx + 1)] - all_dates[Int(midmonth_idx)]).value
         map(
             x ->
                 bcf_info.monthly_fields[x] .= scaling_function(
@@ -192,3 +197,43 @@ function update_midmonth_data!(date, bcf_info)
 end
 
 next_date_in_file(bcfile_info) = bcfile_info.all_dates[bcfile_info.segment_idx[1] + Int(1)]
+
+# IO - daily
+"""
+    interpolate_midmonth_to_daily(date, bcf_info)
+
+Interpolates linearly between two `Fields` in the `bcf_info` struct, or returns the first Field if interpolation is switched off. 
+"""
+function interpolate_midmonth_to_daily(date, bcf_info)
+
+    if bcf_info.interpolate_daily && bcf_info.segment_length[1] > FT(0)
+        segment_length = bcf_info.segment_length
+        segment_idx = bcf_info.segment_idx
+        all_dates = bcf_info.all_dates
+        monthly_fields = bcf_info.monthly_fields
+
+        return interpol.(
+            monthly_fields[1],
+            monthly_fields[2],
+            FT((date - all_dates[Int(segment_idx[1])]).value),
+            FT(segment_length[1]),
+        )
+    else
+        return bcf_info.monthly_fields[1]
+    end
+end
+
+"""
+    interpol(f1::FT, f2::FT, Δt_tt1::FT, Δt_t2t1::FT) where {FT}
+
+Performs linear interpolation of `f` at time `t` within a segment `Δt_t2t1 = (t2 - t1)`, of fields `f1` and `f2`, with `t2 > t1`. 
+
+`Δt_tt1 = (t - t1)`
+`f(t1) = f1`
+
+"""
+function interpol(f1::FT, f2::FT, Δt_tt1::FT, Δt_t2t1::FT) where {FT}
+    interp_fraction = Δt_tt1 / Δt_t2t1
+    @assert abs(interp_fraction) <= FT(1) "time interpolation weights must be <= 1, but `interp_fraction` = $interp_fraction"
+    return f1 * interp_fraction + f2 * (FT(1) - interp_fraction)
+end

experiments/AMIP/moist_mpi_earth/coupler_utils/calendar_timer.jl

@@ -33,7 +33,7 @@ Evaluate `ex` when `model_date` is on/after `callback_date` and do nothing other
 """
 macro calendar_callback(ex::Expr, model_date::Union{Symbol, Expr}, callback_date::Union{Symbol, Expr})
     quote
-        if Dates.days($model_date - $callback_date) < FT(0)
+        if ($model_date - $callback_date).value < FT(0)
             nothing
         else
             eval($ex)

experiments/AMIP/moist_mpi_earth/coupler_utils/unit_tester.jl

@@ -1,7 +1,5 @@
 # unit_tester - contains temporaty unit tests that will be moved to the `/test` folder
 
-# unit test for update_midmonth_data! (move during interface revamp)
-
 import Test: @test
 using Pkg
 using Dates
@@ -16,19 +14,22 @@ using Dates
 # include("general_helper.jl")
 # include("bcfile_reader.jl")
 
-# 1. update_midmonth_data! loop
+# test setup
 FT = Float32
 
 date = date0 = DateTime(1979, 01, 01)
 tspan = (1, 90 * 86400) # Jan-Mar
-Δt = 1 * 86400
+Δt = 1 * 3600
 
 domain = Domains.SphereDomain(FT(6731e3))
 mesh = Meshes.EquiangularCubedSphere(domain, 4)
 topology = Topologies.Topology2D(mesh)
 quad = Spaces.Quadratures.GLL{5}()
 boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad)
 
+# IO monthly
+# unit test for update_midmonth_data!
+
 land_mask_t = ClimaCore.Fields.zeros(boundary_space_t)
 
 _info = bcfile_info_init(
@@ -46,7 +47,6 @@ cs_t = (; _info = _info, date = [date], SST_all = [], updating_dates = [])
 @show "Starting coupling loop"
 function test_update_midmonth_callback()
     # step in time
-    ct = 0
     walltime = @elapsed for t in ((tspan[1] + Δt):Δt:tspan[end])
         cs_t.date[1] = current_date(t) # if not global, `date`` is not updated. Check that this still runs when distributed.
 
@@ -64,3 +64,13 @@ function test_update_midmonth_callback()
 end
 
 test_update_midmonth_callback()
+
+# IO daily
+function test_interpol(boundary_space)
+    f1 = zeros(boundary_space)
+    f2 = ones(boundary_space)
+    out = interpol.(f1, f2, FT(15), FT(30))
+    @test parent(out)[1] ≈ FT(0.5)
+end
+
+test_interpol(boundary_space_t)