From c123e61944bfad81c9fd31e386ef8fc073fda7c7 Mon Sep 17 00:00:00 2001 From: Julia Sloan Date: Tue, 25 Oct 2022 19:50:48 -0700 Subject: [PATCH] BCReader module implemented, rebased --- .buildkite/pipeline.yml | 11 + docs/make.jl | 2 +- docs/src/bcreader.md | 23 ++ src/BCReader.jl | 309 +++++++++++++++++++++++++ src/CallbackManager.jl | 43 ++++ src/ClimaCoupler.jl | 3 + src/Regridder.jl | 6 +- src/Utilities.jl | 4 +- test/TestHelper.jl | 2 +- test/bcreader_tests.jl | 313 ++++++++++++++++++++++++++ test/mpi_tests/bcreader_mpi_tests.jl | 231 +++++++++++++++++++ test/mpi_tests/regridder_mpi_tests.jl | 49 ++-- test/mpi_tests/run_mpi_tests.jl | 1 + test/regridder_tests.jl | 295 ++++++++++++------------ test/runtests.jl | 1 + 15 files changed, 1123 insertions(+), 170 deletions(-) create mode 100644 docs/src/bcreader.md create mode 100644 src/BCReader.jl create mode 100644 src/CallbackManager.jl create mode 100644 test/bcreader_tests.jl create mode 100644 test/mpi_tests/bcreader_mpi_tests.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 683b3aa2ed..7a22123e53 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -77,6 +77,17 @@ steps: slurm_nodes: 3 slurm_tasks_per_node: 1 + - label: "MPI BCReader unit tests" + key: "bcreader_mpi_tests" + command: "mpiexec julia --color=yes --project=test/ test/mpi_tests/bcreader_mpi_tests.jl" + timeout_in_minutes: 20 + env: + CLIMACORE_DISTRIBUTED: "MPI" + agents: + config: cpu + queue: central + slurm_nodes: 3 + slurm_tasks_per_node: 1 - group: "Integration Tests" steps: diff --git a/docs/make.jl b/docs/make.jl index d473722fd5..4bad849e6e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -57,7 +57,7 @@ experiment_pages = [ "AMIP" => map(s -> "generated/amip/$(s)", readdir(joinpath(@__DIR__, "src/generated/amip"))), ] interface_pages = - ["couplerstate.md", "timestepping.md", "regridder.md", "conservation.md", "utilities.md", "testhelper.md"] + ["couplerstate.md", "timestepping.md", "regridder.md", "conservation.md", "utilities.md", "bcreader.md", "testhelper.md"] pages = Any["Home" => "index.md", "Examples" => experiment_pages, "Coupler Interface" => interface_pages] diff --git a/docs/src/bcreader.md b/docs/src/bcreader.md new file mode 100644 index 0000000000..e928d9d50a --- /dev/null +++ b/docs/src/bcreader.md @@ -0,0 +1,23 @@ +# BCReader + +This module coordinates reading of boundary conditions from NetCDF files, +as well as regridding calls and temporal interpolations from +monthly to daily intervals. + +## BCReader API + +```@docs + ClimaCoupler.BCReader.BCFileInfo + ClimaCoupler.BCReader.bcfile_info_init + ClimaCoupler.BCReader.update_midmonth_data! + ClimaCoupler.BCReader.next_date_in_file + ClimaCoupler.BCReader.interpolate_midmonth_to_daily +``` + + +## BCReader Internal Functions + +```@docs + ClimaCoupler.BCReader.no_scaling + ClimaCoupler.BCReader.interpol +``` \ No newline at end of file diff --git a/src/BCReader.jl b/src/BCReader.jl new file mode 100644 index 0000000000..43ee7d64e2 --- /dev/null +++ b/src/BCReader.jl @@ -0,0 +1,309 @@ +#= + BCReader + +This module coordinates reading of boundary conditions from NetCDF files, +as well as regridding calls and temporal interpolations from +monthly to daily intervals. +=# +module BCReader + +using ..Utilities, ..Regridder, ..CallbackManager +using ClimaCore: Fields +using ClimaComms +using Dates +using JLD2 + +export BCFileInfo, bcfile_info_init, update_midmonth_data!, next_date_in_file, interpolate_midmonth_to_daily + + +""" + BCFileInfo + +Stores information specific to each boundary condition from a file and each variable. + +# Inputs + FT::F # float type + bcfile_dir::b # directory of the BC file + comms_ctx::X # communication context used for MPI + hd_outfile_root::S # filename root for regridded data + varname::V # name of the variable + all_dates::D # vector of all dates contained in the original data file + segment_idx::Vector{Int} # index of the monthly data in the file + segment_idx0::Vector{Int} # `segment_idx` of the file data that is closest to date0 + 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::Bool # switch to trigger daily interpolation + land_mask::M # mask with 1 = land, 0 = ocean / sea-ice +""" +struct BCFileInfo{F, B, X, S, V, D, C, O, M} + FT::F + bcfile_dir::B + comms_ctx::X + hd_outfile_root::S + varname::V + all_dates::D + segment_idx::Vector{Int} + segment_idx0::Vector{Int} + monthly_fields::C + segment_length::Vector{Int} + scaling_function::O + interpolate_daily::Bool + land_mask::M +end + +""" + no_scaling(field, bcf_info) + +Remap the values of a field onto the space of the +BCFileInfo's land_mask without scaling. + +# Arguments +- `field`: [Fields.Field] contains the values to be remapped. +- `bcf_info`: [BCFileInfo] contains a land_mask to remap onto the space of. +""" +no_scaling(field::Fields.Field, bcf_info::BCFileInfo) = Utilities.swap_space!(field, axes(bcf_info.land_mask)) + +""" + bcfile_info_init( + FT, + datafile_rll, + varname, + boundary_space; + comms_ctx = ClimaComms.SingletonCommsContext(), + interpolate_daily = false, + segment_idx0 = nothing, + scaling_function = no_scaling, + land_mask = nothing, + date0 = nothing, + mono = true, + ) + +Regrids from lat-lon grid to cgll grid, saving the output in a new file, +and returns the info packaged in a single struct. + +# Arguments +- `FT`: [DataType] Float type. +- `bcfile_dir`: [String] directory the BC file is stored in. +- `datafile_rll`: [String] file containing data to regrid. +- `varname`: [String] name of the variable to be regridded. +- `boundary_space`: [Spaces.AbstractSpace] the space to which we are mapping. +- `comms_ctx`: [ClimaComms.AbstractCommsContext] context used for this operation. +- `interpolate_daily`: [Bool] switch to trigger daily interpolation. +- `segment_idx0`: [Vector{Int}] index of the file data that is closest to date0. +- `scaling function`: [Function] scales, offsets or transforms `varname`. +- `land_mask`: [Fields.field] mask with 1 = land, 0 = ocean / sea-ice. +- `date0`: [Dates.DateTime] start date of the file data. +- `mono`: [Bool] flag for monotone remapping of `datafile_rll`. + +# Returns +- BCFileInfo +""" +function bcfile_info_init( + FT, + bcfile_dir, + datafile_rll, + varname, + boundary_space; + comms_ctx = ClimaComms.SingletonCommsContext(), + interpolate_daily = false, + segment_idx0 = nothing, + scaling_function = no_scaling, + land_mask = nothing, + date0 = nothing, + mono = true, +) + + # regrid all times and save to hdf5 files + hd_outfile_root = varname * "_cgll" + if ClimaComms.iamroot(comms_ctx) + Regridder.hdwrite_regridfile_rll_to_cgll( + FT, + bcfile_dir, + datafile_rll, + varname, + boundary_space; + hd_outfile_root = hd_outfile_root, + mono = mono, + ) + end + ClimaComms.barrier(comms_ctx) + data_dates = JLD2.load(joinpath(bcfile_dir, hd_outfile_root * "_times.jld2"), "times") + + # init time tracking info + current_fields = Fields.zeros(FT, boundary_space), Fields.zeros(FT, boundary_space) + segment_length = [Int(0)] + + # unless the start file date is specified, find the closest one to the start date + segment_idx0 = + segment_idx0 != nothing ? segment_idx0 : + [ + argmin( + abs.( + parse(FT, CallbackManager.datetime_to_strdate(date0)) .- + parse.(FT, CallbackManager.datetime_to_strdate.(data_dates[:])) + ), + ), + ] + + return BCFileInfo( + FT, + bcfile_dir, + comms_ctx, + hd_outfile_root, + varname, + data_dates, + deepcopy(segment_idx0), + segment_idx0, + current_fields, + segment_length, + scaling_function, + interpolate_daily, + land_mask, + ) +end + +# IO - monthly +""" + update_midmonth_data!(date, bcf_info) + +Extracts boundary condition data from regridded (to model grid) NetCDF files. +The times for which data is extracted depends on the specifications in the +`bcf_info` struct). + +# Arguments +- `date`: [Dates.DateTime] start date for data. +- `bcf_info`: [BCFileInfo] containing boundary condition data. +""" +function update_midmonth_data!(date, bcf_info::BCFileInfo) + # monthly count + (; FT, bcfile_dir, comms_ctx, hd_outfile_root, varname, all_dates, scaling_function) = bcf_info + midmonth_idx = bcf_info.segment_idx[1] + midmonth_idx0 = bcf_info.segment_idx0[1] + + if (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])" + bcf_info.monthly_fields[1] .= scaling_function( + Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx0)], varname, comms_ctx), + bcf_info, + ) + bcf_info.monthly_fields[2] .= deepcopy(bcf_info.monthly_fields[1]) + bcf_info.segment_length .= Int(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])" + bcf_info.monthly_fields[1] .= scaling_function( + Regridder.read_from_hdf5( + bcfile_dir, + hd_outfile_root, + all_dates[Int(length(all_dates))], + varname, + comms_ctx, + ), + bcf_info, + ) + bcf_info.monthly_fields[2] .= deepcopy(bcf_info.monthly_fields[1]) + bcf_info.segment_length .= Int(0) + + # throw error when there are closer initial indices for the bc file data that matches this date0 + elseif Dates.days(date - all_dates[Int(midmonth_idx + 1)]) > 2 + nearest_idx = argmin( + abs.( + parse(FT, CallbackManager.datetime_to_strdate(date)) .- + parse.(FT, CallbackManager.datetime_to_strdate.(all_dates[:])) + ), + ) + # TODO test this + bcf_info.segment_idx[1] = midmonth_idx = midmonth_idx0 = nearest_idx + @warn "init data does not correspond to start date. Initializing with `SIC_info.segment_idx[1] = midmonth_idx = midmonth_idx0 = $nearest_idx` for this start date" + + # date crosses to the next month + elseif Dates.days(date - all_dates[Int(midmonth_idx)]) > 0 + 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 + bcf_info.monthly_fields[1] .= scaling_function( + Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx)], varname, comms_ctx), + bcf_info, + ) + bcf_info.monthly_fields[2] .= scaling_function( + Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx + 1)], varname, comms_ctx), + bcf_info, + ) + + else + throw(ErrorException("Check boundary file specification")) + end +end + +""" + next_date_in_file(bcf_info) + +Returns the next date stored in the file `bcfile_info` struct after the +current date index given by `segment_idx`. +Note: this function does not update `segment_idx`, so repeated calls will +return the same value unless `segment_idx` is modified elsewhere in between. + +# Arguments +- `bcf_info`: [BCFileInfo] containing the date information. + +# Returns +- Dates.DateTime +""" +next_date_in_file(bcf_info::BCFileInfo) = bcf_info.all_dates[bcf_info.segment_idx[1] + Int(1)] + +# IO - daily +""" + interpolate_midmonth_to_daily(FT, date, bcf_info::BCFileInfo) + +Interpolates linearly between two `Fields` in the `bcf_info` struct, +or returns the first Field if interpolation is switched off. + +# Arguments +- `FT`: [DataType] Float type. +- `date`: [Dates.DateTime] start date for data. +- `bcf_info`: [BCFileInfo] contains fields to be interpolated. + +# Returns +- Fields.field +""" +function interpolate_midmonth_to_daily(FT, date, bcf_info::BCFileInfo) + if bcf_info.interpolate_daily && bcf_info.segment_length[1] > FT(0) + (; segment_length, segment_idx, all_dates, monthly_fields) = bcf_info + + 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) + +Performs linear interpolation of `f` at time `t` within +a segment `Δt_t2t1 = (t2 - t1)`, of fields `f1` and `f2`, with `t2 > t1`. + +# Arguments +- `f1`: [FT] first value to be interpolated (`f(t1) = f1`). +- `f2`: [FT] second value to be interpolated. +- `Δt_tt1`: [FT] time between `t1` and some `t` (`Δt_tt1 = (t - t1)`). +- `Δt_t2t1`: [FT] time between `t1` and `t2`. + +# Returns +- FT +""" +function interpol(f1::FT, f2::FT, Δt_tt1::FT, Δt_t2t1::FT) where {FT} + @assert Δt_t2t1 > FT(0) "t2 must be > t1, but `Δt_t2t1` = $Δt_t2t1" + 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 + +end diff --git a/src/CallbackManager.jl b/src/CallbackManager.jl new file mode 100644 index 0000000000..7777927937 --- /dev/null +++ b/src/CallbackManager.jl @@ -0,0 +1,43 @@ +#= + CallbackManager - TODO WIP + +This module facilitates calendar functions and temporal interpolations +of data. +=# + +module CallbackManager + +using ..Utilities +using Dates + +export current_date, strdate_to_datetime, datetime_to_strdate, calendar_callback + +# TODO remove calls to calendar_callback - can be replaced with `<` + +""" + current_date(cs, t) + +Return the model date +""" +current_date(cs, t) = cs.dates.date0[1] + Dates.Second(t) + +""" + strdate_to_datetime(strdate) + +Convert from String ("YYYYMMDD") to Date format +""" +strdate_to_datetime(strdate::String) = + Dates.DateTime(parse(Int, strdate[1:4]), parse(Int, strdate[5:6]), parse(Int, strdate[7:8])) # required by the official AMIP input files + +""" + datetime_to_strdate(datetime) + +Convert from Date to String ("YYYYMMDD") format +""" +datetime_to_strdate(datetime::DateTime) = + string(Dates.year(datetime)) * + string(string(lpad(Dates.month(datetime), 2, "0"))) * + string(lpad(Dates.day(datetime), 2, "0")) + + +end diff --git a/src/ClimaCoupler.jl b/src/ClimaCoupler.jl index cc5bb246c3..6be3be8451 100644 --- a/src/ClimaCoupler.jl +++ b/src/ClimaCoupler.jl @@ -12,5 +12,8 @@ include("../test/TestHelper.jl") include("Utilities.jl") include("Regridder.jl") include("ConservationChecker.jl") +include("CallbackManager.jl") +include("BCReader.jl") + end diff --git a/src/Regridder.jl b/src/Regridder.jl index 84b2ad0ddf..2a3ec7a399 100644 --- a/src/Regridder.jl +++ b/src/Regridder.jl @@ -8,7 +8,7 @@ via ClimaCoreTempestRemap wrappers. module Regridder using ClimaCoupler.Utilities -using ClimaCore: ClimaCore, Meshes, Domains, Topologies, Spaces, Fields, InputOutput +using ClimaCore: Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaComms using NCDatasets using ClimaCoreTempestRemap @@ -102,7 +102,7 @@ function hdwrite_regridfile_rll_to_cgll( topology = Topologies.Topology2D(space.topology.mesh, Topologies.spacefillingcurve(space.topology.mesh)) Nq = Spaces.Quadratures.polynomial_degree(space.quadrature_style) + 1 - space_undistributed = ClimaCore.Spaces.SpectralElementSpace2D(topology, ClimaCore.Spaces.Quadratures.GLL{Nq}()) + space_undistributed = Spaces.SpectralElementSpace2D(topology, Spaces.Quadratures.GLL{Nq}()) if isfile(datafile_cgll) == false isdir(REGRID_DIR) ? nothing : mkpath(REGRID_DIR) @@ -375,7 +375,7 @@ function land_sea_mask( ) end ClimaComms.barrier(comms_ctx) - file_dates = load(joinpath(REGRID_DIR, outfile_root * "_times.jld2"), "times") + file_dates = JLD2.load(joinpath(REGRID_DIR, outfile_root * "_times.jld2"), "times") mask = read_from_hdf5(REGRID_DIR, outfile_root, file_dates[1], varname, comms_ctx) mask = swap_space!(mask, boundary_space) # needed if we are reading from previous run return mono ? mask : binary_mask.(mask, threshold = threshold) diff --git a/src/Utilities.jl b/src/Utilities.jl index 133cf1c196..94999b6fe1 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -6,11 +6,11 @@ modules in the coupler. =# module Utilities -using ClimaCore: ClimaCore, Fields, Spaces, Domains, Meshes, Topologies -using ClimaComms +using ClimaCore: Fields, Spaces export CoupledSimulation, heaviside, swap_space!, create_space + """ Stores information needed to run a simulation with the coupler. """ diff --git a/test/TestHelper.jl b/test/TestHelper.jl index f057c2d9d4..a75cc20097 100644 --- a/test/TestHelper.jl +++ b/test/TestHelper.jl @@ -6,7 +6,7 @@ various files in the test folder. =# module TestHelper -using ClimaCore: ClimaCore, Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput +using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaComms using NCDatasets diff --git a/test/bcreader_tests.jl b/test/bcreader_tests.jl new file mode 100644 index 0000000000..a97c412d08 --- /dev/null +++ b/test/bcreader_tests.jl @@ -0,0 +1,313 @@ +#= + Unit tests for ClimaCoupler BCReader module +=# + +using ClimaCoupler: Regridder, BCReader, CallbackManager +using ClimaCore: Fields, Meshes, Domains, Topologies, Spaces +using ClimaComms +using Test +using Dates +using NCDatasets +import ArtifactWrappers as AW + +# get the paths to the necessary data files - sst map, land sea mask +include(joinpath(@__DIR__, "..", "artifacts", "artifact_funcs.jl")) +sst_data = joinpath(sst_dataset_path(), "sst.nc") +mask_data = joinpath(mask_dataset_path(), "seamask.nc") + +for FT in (Float32, Float64) + @testset "test next_date_in_file for FT=$FT" begin + dummy_dates = Vector(range(DateTime(1999, 1, 1); step = Day(1), length = 10)) + date0 = dummy_dates[1] + segment_idx0 = [ + argmin( + abs.( + parse(FT, CallbackManager.datetime_to_strdate(date0)) .- + parse.(FT, CallbackManager.datetime_to_strdate.(dummy_dates[:])) + ), + ), + ] + + bcf_info = BCReader.BCFileInfo( + FT, # FT + "", # bcfile_dir + ClimaComms.SingletonCommsContext(), # comms_ctx + "", # hd_outfile_root + "", # varname + dummy_dates, # all_dates + deepcopy(segment_idx0), # segment_idx + segment_idx0, # segment_idx0 + nothing, # monthly_fields + Int[], # segment_length + nothing, # scaling_function + false, # interpolate_daily + nothing, # land_mask + ) + + idx = segment_idx0[1] + current_date = date0 + next_date = BCReader.next_date_in_file(bcf_info) + @test current_date == dummy_dates[idx] + + for i in 1:(length(dummy_dates) - 2) + current_date = next_date + bcf_info.segment_idx[1] += Int(1) + next_date = BCReader.next_date_in_file(bcf_info) + idx = segment_idx0[1] + i + @test next_date == dummy_dates[idx + 1] + end + end + + @testset "test interpolate_midmonth_to_daily for FT=$FT" begin + # test interpolate_midmonth_to_daily with interpolation + interpolate_daily = true + dummy_dates = Vector(range(DateTime(1999, 1, 1); step = Day(1), length = 100)) + segment_idx0 = [Int(1)] + + # these values give an `interp_fraction` of 0.5 in `interpol` for ease of testing + date0 = dummy_dates[Int(segment_idx0[1] + 1)] + segment_length = [Int(2) * ((date0 - dummy_dates[Int(segment_idx0[1])]).value)] + + radius = FT(6731e3) + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, 4) + topology = Topologies.Topology2D(mesh) + quad = Spaces.Quadratures.GLL{Nq}() + boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad) + monthly_fields = (zeros(boundary_space_t), ones(boundary_space_t)) + + bcf_info_interp = BCReader.BCFileInfo( + FT, # FT + "", # bcfile_dir + ClimaComms.SingletonCommsContext(), # comms_ctx + "", # hd_outfile_root + "", # varname + dummy_dates, # all_dates + deepcopy(segment_idx0), # segment_idx + segment_idx0, # segment_idx0 + monthly_fields, # monthly_fields + segment_length, # segment_length + nothing, # scaling_function + interpolate_daily, # interpolate_daily + nothing, # land_mask + ) + @test BCReader.interpolate_midmonth_to_daily(FT, date0, bcf_info_interp) == ones(boundary_space_t) .* FT(0.5) + + # test interpolate_midmonth_to_daily without interpolation + interpolate_daily = false + + bcf_info_no_interp = BCReader.BCFileInfo( + FT, # FT + "", # bcfile_dir + ClimaComms.SingletonCommsContext(), # comms_ctx + "", # hd_outfile_root + "", # varname + dummy_dates, # all_dates + deepcopy(segment_idx0), # segment_idx + segment_idx0, # segment_idx0 + monthly_fields, # monthly_fields + segment_length, # segment_length + nothing, # scaling_function + interpolate_daily, # interpolate_daily + nothing, # land_mask + ) + @test BCReader.interpolate_midmonth_to_daily(FT, date0, bcf_info_no_interp) == monthly_fields[1] + end + + # Add tests which use TempestRemap here - + # TempestRemap is not built on Windows because of NetCDF support limitations + # `bcf_info_init` uses TR via a call to `hdwrite_regridfile_rll_to_cgll` + if !Sys.iswindows() + @testset "test update_midmonth_data! for FT=$FT" begin + # setup for test + date0 = date1 = DateTime(1979, 01, 01, 01, 00, 00) + date = DateTime(1979, 01, 01, 00, 00, 00) + tspan = (1, 90 * 86400) # Jan-Mar + Δt = 1 * 3600 + + radius = FT(6731e3) + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, 4) + topology = Topologies.Topology2D(mesh) + quad = Spaces.Quadratures.GLL{Nq}() + boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad) + + land_mask_t = Fields.zeros(boundary_space_t) + dummy_data = (; test_data = zeros(axes(land_mask_t))) + + datafile_rll = sst_data + varname = "SST" + + regrid_dir = "bcreader_regrid_dir" + isdir(regrid_dir) ? nothing : mkpath(regrid_dir) + + bcf_info = BCReader.bcfile_info_init( + FT, + regrid_dir, + datafile_rll, + varname, + boundary_space_t, + segment_idx0 = [Int(1309)], + interpolate_daily = false, + land_mask = land_mask_t, + ) + + dates = (; date = [date], date0 = [date0], date1 = [date1]) + cs_t = (; + _info = bcf_info, + dates = dates, + SST_all = [], + updating_dates = [], + monthly_state_diags = (; fields = deepcopy(dummy_data), ct = [0]), + ) + + # step in time + walltime = @elapsed for t in ((tspan[1] + Δt):Δt:tspan[end]) + cs_t.dates.date[1] = CallbackManager.current_date(cs_t, t) # if not global, `date`` is not updated. Check that this still runs when distributed. + + model_date = cs_t.dates.date[1] + callback_date = BCReader.next_date_in_file(cs_t._info) + + # TODO investigate if macro would be faster here + if (model_date >= callback_date) + BCReader.update_midmonth_data!(model_date, cs_t._info) + push!(cs_t.SST_all, deepcopy(cs_t._info.monthly_fields[1])) + push!(cs_t.updating_dates, deepcopy(model_date)) + end + + end + + # test if the SST field was modified + @test cs_t.SST_all[end] !== cs_t.SST_all[end - 1] + # check that the final file date is as expected + @test Date(cs_t.updating_dates[end]) == Date(1979, 03, 16) + + # test warning/error cases + current_fields = Fields.zeros(FT, boundary_space_t), Fields.zeros(FT, boundary_space_t) + + # use this function to reset values between test cases + function reset_bcf_info(bcf_info) + bcf_info.monthly_fields[1] .= current_fields[1] + bcf_info.monthly_fields[2] .= current_fields[2] + bcf_info.segment_length[1] = Int(1) + end + + hd_outfile_root = varname * "_cgll" + comms_ctx = ClimaComms.SingletonCommsContext() + + # case 1: segment_idx == segment_idx0, date < all_dates[segment_idx] + bcf_info.segment_idx[1] = bcf_info.segment_idx0[1] + date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1)) + BCReader.update_midmonth_data!(date, bcf_info) + + @test bcf_info.monthly_fields[1] == bcf_info.scaling_function( + Regridder.read_from_hdf5( + regrid_dir, + hd_outfile_root, + bcf_info.all_dates[Int(bcf_info.segment_idx0[1])], + varname, + comms_ctx, + ), + bcf_info, + ) + @test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1] + @test bcf_info.segment_length[1] == Int(0) + + # case 2: date > all_dates[end - 1] + reset_bcf_info(bcf_info) + date = DateTime(bcf_info.all_dates[end - 1] + Dates.Day(1)) + BCReader.update_midmonth_data!(date, bcf_info) + + @test bcf_info.monthly_fields[1] == bcf_info.scaling_function( + Regridder.read_from_hdf5( + regrid_dir, + hd_outfile_root, + bcf_info.all_dates[Int(length(bcf_info.all_dates))], + varname, + comms_ctx, + ), + bcf_info, + ) + @test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1] + @test bcf_info.segment_length[1] == Int(0) + + # case 3: date - all_dates[segment_idx + 1] > 2 + reset_bcf_info(bcf_info) + date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1] + 1] + Dates.Day(3)) + BCReader.update_midmonth_data!(date, bcf_info) + + nearest_idx = argmin( + abs.( + parse(FT, CallbackManager.datetime_to_strdate(date)) .- + parse.(FT, CallbackManager.datetime_to_strdate.(bcf_info.all_dates[:])) + ), + ) + @test bcf_info.segment_idx[1] == bcf_info.segment_idx0[1] == nearest_idx + + # case 4: everything else + reset_bcf_info(bcf_info) + bcf_info.segment_idx[1] = bcf_info.segment_idx0[1] + Int(1) + date = bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1) + + @test_throws ErrorException BCReader.update_midmonth_data!(date, bcf_info) + + rm(regrid_dir; recursive = true, force = true) + end + + @testset "test bcf_info_init for FT=$FT" begin + # setup for test + radius = FT(6731e3) + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, 4) + topology = Topologies.Topology2D(mesh) + quad = Spaces.Quadratures.GLL{Nq}() + boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad) + land_mask_t = Fields.zeros(boundary_space_t) + + datafile_rll = mask_data + varname = "LSMASK" + mono = true + + regrid_dir = "bcreader_regrid_dir" + isdir(regrid_dir) ? nothing : mkpath(regrid_dir) + + bcf_info = BCReader.bcfile_info_init( + FT, + regrid_dir, + datafile_rll, + varname, + boundary_space_t, + segment_idx0 = [Int(1309)], + land_mask = land_mask_t, + mono = mono, + ) + + # test that created object exists and has correct components + @test @isdefined(bcf_info) + @test all(parent(bcf_info.land_mask) .== 0) + + # construct weightfile name to test values + hd_outfile_root = varname * "_cgll" + outfile = hd_outfile_root * ".nc" + outfile_root = mono ? outfile[1:(end - 3)] * "_mono" : outfile[1:(end - 3)] + weightfile = joinpath(regrid_dir, outfile_root * "_remap_weights.nc") + + # test monotone remapping (all weights in [0, 1]) + nt = NCDataset(weightfile) do weights + max_weight = maximum(weights["S"]) + min_weight = minimum(weights["S"]) + (; max_weight, min_weight) + end + (; max_weight, min_weight) = nt + + @test max_weight <= FT(1.0) || isapprox(max_weight, FT(1.0), atol = 1e-16) + @test min_weight >= FT(0.0) || isapprox(min_weight, FT(0.0), atol = 1e-16) + + # delete testing directory and files + rm(regrid_dir; recursive = true, force = true) + end + end +end diff --git a/test/mpi_tests/bcreader_mpi_tests.jl b/test/mpi_tests/bcreader_mpi_tests.jl new file mode 100644 index 0000000000..ac5dd906d6 --- /dev/null +++ b/test/mpi_tests/bcreader_mpi_tests.jl @@ -0,0 +1,231 @@ +#= + Unit tests for ClimaCoupler BCReader module functions which require MPI + +These are in a separate testing file from the other BCReader unit tests so +that MPI can be enabled for testing of these functions. +=# + +using ClimaCoupler: BCReader, CallbackManager, Regridder +using ClimaCore: Fields, Meshes, Domains, Topologies, Spaces +using ClimaComms +using Test +using Dates +using NCDatasets +import ArtifactWrappers as AW +using ClimaCommsMPI + +# Get the path to the necessary data file - sst map +include(joinpath(@__DIR__, "..", "..", "artifacts", "artifact_funcs.jl")) +sst_data = joinpath(sst_dataset_path(), "sst.nc") + +# set up MPI communications context +comms_ctx = ClimaCommsMPI.MPICommsContext() +pid, nprocs = ClimaComms.init(comms_ctx) +ClimaComms.barrier(comms_ctx) + +@testset "test bcf_info_init with MPI" begin + for FT in (Float32, Float64) + # setup for test + radius = FT(6731e3) + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, 4) + topology = Topologies.DistributedTopology2D(comms_ctx, mesh, Topologies.spacefillingcurve(mesh)) + quad = Spaces.Quadratures.GLL{Nq}() + boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad) + land_mask_t = Fields.zeros(boundary_space_t) + + datafile_rll = sst_data + varname = "SST" + mono = true + + regrid_dir = "bcreader_regrid_dir" + isdir(regrid_dir) ? nothing : mkpath(regrid_dir) + + bcf_info = BCReader.bcfile_info_init( + FT, + regrid_dir, + datafile_rll, + varname, + boundary_space_t, + comms_ctx = comms_ctx, + segment_idx0 = [Int(1309)], + land_mask = land_mask_t, + mono = mono, + ) + + # test that created object exists and has correct components + @test @isdefined(bcf_info) + @test all(parent(bcf_info.land_mask) .== 0) + + # construct weightfile name to test values + hd_outfile_root = varname * "_cgll" + outfile = hd_outfile_root * ".nc" + outfile_root = mono ? outfile[1:(end - 3)] * "_mono" : outfile[1:(end - 3)] + weightfile = joinpath(regrid_dir, outfile_root * "_remap_weights.nc") + + # test monotone remapping (all weights in [0, 1]) + nt = NCDataset(weightfile) do weights + max_weight = maximum(weights["S"]) + min_weight = minimum(weights["S"]) + (; max_weight, min_weight) + end + (; max_weight, min_weight) = nt + + @test max_weight <= FT(1.0) || isapprox(max_weight, FT(1.0), atol = 1e-16) + @test min_weight >= FT(0.0) || isapprox(min_weight, FT(0.0), atol = 1e-16) + + # delete testing directory and files + rm(regrid_dir; recursive = true, force = true) + end +end + +@testset "test update_midmonth_data! with MPI" begin + for FT in (Float32, Float64) + # setup for test + date0 = date1 = DateTime(1979, 01, 01, 01, 00, 00) + date = DateTime(1979, 01, 01, 00, 00, 00) + tspan = (1, 90 * 86400) # Jan-Mar + Δt = 1 * 3600 + + radius = FT(6731e3) + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, 4) + topology = Topologies.DistributedTopology2D(comms_ctx, mesh, Topologies.spacefillingcurve(mesh)) + quad = Spaces.Quadratures.GLL{Nq}() + boundary_space_t = Spaces.SpectralElementSpace2D(topology, quad) + + land_mask_t = Fields.zeros(boundary_space_t) + dummy_data = (; test_data = zeros(axes(land_mask_t))) + + datafile_rll = sst_data + varname = "SST" + + regrid_dir = "bcreader_regrid_dir" + isdir(regrid_dir) ? nothing : mkpath(regrid_dir) + + ClimaComms.barrier(comms_ctx) + + bcf_info = BCReader.bcfile_info_init( + FT, + regrid_dir, + datafile_rll, + varname, + boundary_space_t, + comms_ctx = comms_ctx, + segment_idx0 = [Int(1309)], + interpolate_daily = false, + land_mask = land_mask_t, + ) + + dates = (; date = [date], date0 = [date0], date1 = [date1]) + cs_t = (; + _info = bcf_info, + dates = dates, + SST_all = [], + updating_dates = [], + monthly_state_diags = (; fields = deepcopy(dummy_data), ct = [0]), + ) + + ClimaComms.barrier(comms_ctx) + + # step in time + walltime = @elapsed for t in ((tspan[1] + Δt):Δt:tspan[end]) + cs_t.dates.date[1] = CallbackManager.current_date(cs_t, t) # if not global, `date`` is not updated. Check that this still runs when distributed. + + model_date = cs_t.dates.date[1] + callback_date = BCReader.next_date_in_file(cs_t._info) + + # TODO investigate if macro would be faster here + if (model_date >= callback_date) + BCReader.update_midmonth_data!(model_date, cs_t._info) + push!(cs_t.SST_all, deepcopy(cs_t._info.monthly_fields[1])) + push!(cs_t.updating_dates, deepcopy(model_date)) + end + + end + + ClimaComms.barrier(comms_ctx) + # test if the SST field was modified + @test cs_t.SST_all[end] !== cs_t.SST_all[end - 1] + # check that the final file date is as expected + @test Date(cs_t.updating_dates[end]) == Date(1979, 03, 16) + + # test warning/error cases + current_fields = Fields.zeros(FT, boundary_space_t), Fields.zeros(FT, boundary_space_t) + + # use this function to reset values between test cases + function reset_bcf_info(bcf_info) + bcf_info.monthly_fields[1] .= current_fields[1] + bcf_info.monthly_fields[2] .= current_fields[2] + bcf_info.segment_length[1] = Int(1) + end + + hd_outfile_root = varname * "_cgll" + + # case 1: segment_idx == segment_idx0, date < all_dates[segment_idx] + bcf_info.segment_idx[1] = bcf_info.segment_idx0[1] + date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1)) + BCReader.update_midmonth_data!(date, bcf_info) + + ClimaComms.barrier(comms_ctx) + @test bcf_info.monthly_fields[1] == bcf_info.scaling_function( + Regridder.read_from_hdf5( + regrid_dir, + hd_outfile_root, + bcf_info.all_dates[Int(bcf_info.segment_idx0[1])], + varname, + comms_ctx, + ), + bcf_info, + ) + @test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1] + @test bcf_info.segment_length[1] == Int(0) + + # case 2: date > all_dates[end - 1] + reset_bcf_info(bcf_info) + date = DateTime(bcf_info.all_dates[end - 1] + Dates.Day(1)) + BCReader.update_midmonth_data!(date, bcf_info) + + ClimaComms.barrier(comms_ctx) + + @test bcf_info.monthly_fields[1] == bcf_info.scaling_function( + Regridder.read_from_hdf5( + regrid_dir, + hd_outfile_root, + bcf_info.all_dates[Int(length(bcf_info.all_dates))], + varname, + comms_ctx, + ), + bcf_info, + ) + @test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1] + @test bcf_info.segment_length[1] == Int(0) + + # case 3: date - all_dates[segment_idx + 1] > 2 + reset_bcf_info(bcf_info) + date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1] + 1] + Dates.Day(3)) + BCReader.update_midmonth_data!(date, bcf_info) + + nearest_idx = argmin( + abs.( + parse(FT, CallbackManager.datetime_to_strdate(date)) .- + parse.(FT, CallbackManager.datetime_to_strdate.(bcf_info.all_dates[:])) + ), + ) + + ClimaComms.barrier(comms_ctx) + @test bcf_info.segment_idx[1] == bcf_info.segment_idx0[1] == nearest_idx + + # case 4: everything else + reset_bcf_info(bcf_info) + bcf_info.segment_idx[1] = bcf_info.segment_idx0[1] + Int(1) + date = bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1) + + ClimaComms.barrier(comms_ctx) + @test_throws ErrorException BCReader.update_midmonth_data!(date, bcf_info) + + rm(regrid_dir; recursive = true, force = true) + end +end diff --git a/test/mpi_tests/regridder_mpi_tests.jl b/test/mpi_tests/regridder_mpi_tests.jl index 493b5042ce..012db1ee74 100644 --- a/test/mpi_tests/regridder_mpi_tests.jl +++ b/test/mpi_tests/regridder_mpi_tests.jl @@ -1,37 +1,40 @@ #= Unit tests for ClimaCoupler Regridder module functions which require MPI -These are in a separate testing file so that MPI can be enabled for -testing of these functions +These are in a separate testing file from the other Regridder unit tests so +that MPI can be enabled for testing of these functions. =# -using ClimaCoupler, ClimaCoupler.TestHelper, ClimaCoupler.Regridder -using ClimaCore: ClimaCore, Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput +using ClimaCoupler: TestHelper, Regridder +using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaCommsMPI using ClimaComms using Dates using Test -FT = Float64 REGRID_DIR = @isdefined(REGRID_DIR) ? REGRID_DIR : joinpath("", "regrid_tmp/") -@testset "test write_to_hdf5 and read_from_hdf5 with MPI" begin - # Set up testing directory - mkpath(REGRID_DIR) - - comms_ctx = ClimaCommsMPI.MPICommsContext() - pid, nprocs = ClimaComms.init(comms_ctx) - - hd_outfile_root = "hdf5_out_test" - tx = Dates.DateTime(1979, 01, 01, 01, 00, 00) - test_space = create_space(FT, comms_ctx = comms_ctx) - input_field = Fields.ones(test_space) - varname = "testdata" +# Set up MPI communications context +# Note that runs will hang if a context is initialized twice in the same file, +# so this context should be shared among all tests in this file. +comms_ctx = ClimaCommsMPI.MPICommsContext() +pid, nprocs = ClimaComms.init(comms_ctx) - @show comms_ctx - ClimaComms.barrier(comms_ctx) - write_to_hdf5(REGRID_DIR, hd_outfile_root, tx, input_field, varname, comms_ctx) - - output_field = read_from_hdf5(REGRID_DIR, hd_outfile_root, tx, varname, comms_ctx) - @test parent(input_field) == parent(output_field) +@testset "test write_to_hdf5 and read_from_hdf5 with MPI" begin + for FT in (Float32, Float64) + # Set up testing directory + mkpath(REGRID_DIR) + + hd_outfile_root = "hdf5_out_test" + tx = Dates.DateTime(1979, 01, 01, 01, 00, 00) + test_space = TestHelper.create_space(FT, comms_ctx = comms_ctx) + input_field = Fields.ones(test_space) + varname = "testdata" + + ClimaComms.barrier(comms_ctx) + Regridder.write_to_hdf5(REGRID_DIR, hd_outfile_root, tx, input_field, varname, comms_ctx) + + output_field = Regridder.read_from_hdf5(REGRID_DIR, hd_outfile_root, tx, varname, comms_ctx) + @test parent(input_field) == parent(output_field) + end end diff --git a/test/mpi_tests/run_mpi_tests.jl b/test/mpi_tests/run_mpi_tests.jl index 7466a951b4..49b92e559b 100644 --- a/test/mpi_tests/run_mpi_tests.jl +++ b/test/mpi_tests/run_mpi_tests.jl @@ -12,4 +12,5 @@ end if !Sys.iswindows() runmpi(joinpath(@__DIR__, "regridder_mpi_tests.jl"), ntasks = 2) + runmpi(joinpath(@__DIR__, "bcreader_mpi_tests.jl"), ntasks = 2) end diff --git a/test/regridder_tests.jl b/test/regridder_tests.jl index eb95791f57..f0001fa44c 100644 --- a/test/regridder_tests.jl +++ b/test/regridder_tests.jl @@ -3,162 +3,177 @@ =# using ClimaCoupler: Utilities, Regridder, TestHelper -using ClimaCore: ClimaCore, Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput +using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaComms using Test using NCDatasets using Dates REGRID_DIR = @isdefined(REGRID_DIR) ? REGRID_DIR : joinpath("", "regrid_tmp/") -FT = Float64 -@testset "test dummmy_remap!" begin - test_space = TestHelper.create_space(FT) - test_field_ones = Fields.ones(test_space) - target_field = Fields.zeros(test_space) - - Regridder.dummmy_remap!(target_field, test_field_ones) - @test parent(target_field) == parent(test_field_ones) -end - -@testset "test update_masks!" begin - test_space = TestHelper.create_space(FT) - # Construct land mask of 0s in top half, 1s in bottom half - land_mask = Fields.ones(test_space) - dims = size(parent(land_mask)) - m = dims[1] - n = dims[2] - parent(land_mask)[1:(m ÷ 2), :, :, :] .= FT(0) - - # Construct ice mask of 0s on left, 0.5s on right - ice_d = Fields.zeros(test_space) - parent(ice_d)[:, (n ÷ 2 + 1):n, :, :] .= FT(0.5) - - # Fill in only the necessary parts of the simulation - cs = Utilities.CoupledSimulation( - nothing, - nothing, - nothing, - nothing, - nothing, - FT, - (; land = land_mask, ice = Fields.zeros(test_space), ocean = Fields.zeros(test_space)), - (;), - (; ice_sim = (; integrator = (; p = (; ice_mask = ice_d)))), - (;), - nothing, - (;), - (;), - (;), - ) - - Regridder.update_masks!(cs) - - # Test that sum of masks is 1 everywhere - @test all(parent(cs.surface_masks.ice .+ cs.surface_masks.land .+ cs.surface_masks.ocean) .== FT(1)) -end - -@testset "test combine_surfaces!" begin - test_space = TestHelper.create_space(FT) - combined_field = Fields.ones(test_space) - - # Initialize weights (masks) and initial values (fields) - masks = (a = 0.0, b = 0.5, c = 2.0, d = -10.0) - fields = (a = 1.0, b = 1.0, c = 1.0, d = 1.0) - - Regridder.combine_surfaces!(combined_field::Fields.Field, masks::NamedTuple, fields::NamedTuple) - @test all(parent(combined_field) .== FT(sum(masks) * sum(fields) / length(fields))) -end - -# Add tests which use TempestRemap here - -# TempestRemap is not built on Windows because of NetCDF support limitations -if !Sys.iswindows() - @testset "test write_to_hdf5 and read_from_hdf5" begin - # Set up testing directory - ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing - mkpath(REGRID_DIR) - - hd_outfile_root = "hdf5_out_test" - tx = Dates.DateTime(1979, 01, 01, 01, 00, 00) +for FT in (Float32, Float64) + @testset "test dummmy_remap!" begin test_space = TestHelper.create_space(FT) - input_field = Fields.ones(test_space) - varname = "testdata" - comms_ctx = ClimaComms.SingletonCommsContext() - - Regridder.write_to_hdf5(REGRID_DIR, hd_outfile_root, tx, input_field, varname, comms_ctx) - - output_field = Regridder.read_from_hdf5(REGRID_DIR, hd_outfile_root, tx, varname, comms_ctx) - @test parent(input_field) == parent(output_field) + test_field_ones = Fields.ones(test_space) + target_field = Fields.zeros(test_space) - # Delete testing directory and files - rm(REGRID_DIR; recursive = true, force = true) + Regridder.dummmy_remap!(target_field, test_field_ones) + @test parent(target_field) == parent(test_field_ones) end - @testset "test remap_field_cgll_to_rll" begin - # Set up testing directory - remap_tmpdir = joinpath(REGRID_DIR, "cgll_to_rll") - ispath(remap_tmpdir) ? rm(remap_tmpdir; recursive = true, force = true) : nothing - mkpath(remap_tmpdir) - name = "testdata" - datafile_rll = remap_tmpdir * "/" * name * "_rll.nc" - + @testset "test update_masks!" begin test_space = TestHelper.create_space(FT) - field = Fields.ones(test_space) + # Construct land mask of 0s in top half, 1s in bottom half + land_mask = Fields.ones(test_space) + dims = size(parent(land_mask)) + m = dims[1] + n = dims[2] + parent(land_mask)[1:(m ÷ 2), :, :, :] .= FT(0) + + # Construct ice mask of 0s on left, 0.5s on right + ice_d = Fields.zeros(test_space) + parent(ice_d)[:, (n ÷ 2 + 1):n, :, :] .= FT(0.5) + + # Fill in only the necessary parts of the simulation + cs = Utilities.CoupledSimulation( + nothing, + nothing, + nothing, + nothing, + nothing, + FT, + (; land = land_mask, ice = Fields.zeros(test_space), ocean = Fields.zeros(test_space)), + (;), + (; ice_sim = (; integrator = (; p = (; ice_mask = ice_d)))), + (;), + nothing, + (;), + (;), + (;), + ) + + Regridder.update_masks!(cs) + + # Test that sum of masks is 1 everywhere + @test all(parent(cs.surface_masks.ice .+ cs.surface_masks.land .+ cs.surface_masks.ocean) .== FT(1)) + end - Regridder.remap_field_cgll_to_rll(name, field, remap_tmpdir, datafile_rll) + @testset "test combine_surfaces!" begin + test_space = TestHelper.create_space(FT) + combined_field = Fields.ones(test_space) - ncdataset_rll = NCDataset(datafile_rll) - @test maximum(ncdataset_rll[name]) == maximum(field) - @test minimum(ncdataset_rll[name]) == minimum(field) + # Initialize weights (masks) and initial values (fields) + masks = (a = 0.0, b = 0.5, c = 2.0, d = -10.0) + fields = (a = 1.0, b = 1.0, c = 1.0, d = 1.0) - # Delete testing directory and files - rm(REGRID_DIR; recursive = true, force = true) + Regridder.combine_surfaces!(combined_field::Fields.Field, masks::NamedTuple, fields::NamedTuple) + @test all(parent(combined_field) .== FT(sum(masks) * sum(fields) / length(fields))) end - @testset "test land_sea_mask" begin - # Test setup - R = FT(6371e3) - test_space = TestHelper.create_space(FT, R = R) - comms_ctx = ClimaComms.SingletonCommsContext() - ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing - mkpath(REGRID_DIR) - - # Initialize dataset of all ones - data_path = joinpath(REGRID_DIR, "ls_mask_data.nc") - varname = "test_data" - TestHelper.gen_ncdata(FT, data_path, varname, FT(1)) - - # Test monotone masking - land_mask_mono = Regridder.land_sea_mask(FT, REGRID_DIR, comms_ctx, data_path, varname, test_space, mono = true) - - # Test no new extrema are introduced in monotone remapping - dataset = NCDataset(data_path) - @test maximum(dataset[varname]) >= maximum(land_mask_mono) && - minimum(dataset[varname]) <= minimum(land_mask_mono) - # Test that monotone remapping a dataset of all ones conserves surface area - @test sum(land_mask_mono) - 4 * π * (R^2) < 10e-14 - - # Delete testing directory and files - rm(REGRID_DIR; recursive = true, force = true) - - # Set up testing directory - ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing - mkpath(REGRID_DIR) - - # Initialize dataset of all 0.5s - data_path = joinpath(REGRID_DIR, "ls_mask_data.nc") - varname = "test_data_halves" - TestHelper.gen_ncdata(FT, data_path, varname, FT(0.5)) - - # Test non-monotone masking - land_mask_halves = - Regridder.land_sea_mask(FT, REGRID_DIR, comms_ctx, data_path, varname, test_space, mono = false) - dataset = NCDataset(data_path) - - # Masking of values below threshold should result in 0 - @test all(parent(land_mask_halves) .== FT(0)) - - # Delete testing directory and files - rm(REGRID_DIR; recursive = true, force = true) + # Add tests which use TempestRemap here - + # TempestRemap is not built on Windows because of NetCDF support limitations + if !Sys.iswindows() + @testset "test write_to_hdf5 and read_from_hdf5" begin + # Set up testing directory + ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing + mkpath(REGRID_DIR) + + hd_outfile_root = "hdf5_out_test" + tx = Dates.DateTime(1979, 01, 01, 01, 00, 00) + test_space = TestHelper.create_space(FT) + input_field = Fields.ones(test_space) + varname = "testdata" + comms_ctx = ClimaComms.SingletonCommsContext() + + Regridder.write_to_hdf5(REGRID_DIR, hd_outfile_root, tx, input_field, varname, comms_ctx) + + output_field = Regridder.read_from_hdf5(REGRID_DIR, hd_outfile_root, tx, varname, comms_ctx) + @test parent(input_field) == parent(output_field) + + # Delete testing directory and files + rm(REGRID_DIR; recursive = true, force = true) + end + + @testset "test remap_field_cgll_to_rll for FT=$FT" begin + # Set up testing directory + remap_tmpdir = joinpath(REGRID_DIR, "cgll_to_rll") + ispath(remap_tmpdir) ? rm(remap_tmpdir; recursive = true, force = true) : nothing + mkpath(remap_tmpdir) + name = "testdata" + datafile_rll = remap_tmpdir * "/" * name * "_rll.nc" + + test_space = TestHelper.create_space(FT) + field = Fields.ones(test_space) + + Regridder.remap_field_cgll_to_rll(name, field, remap_tmpdir, datafile_rll) + + # Test no new extrema are introduced in monotone remapping + nt = NCDataset(datafile_rll) do ds + max_remapped = maximum(ds[name]) + min_remapped = minimum(ds[name]) + (; max_remapped, min_remapped) + end + (; max_remapped, min_remapped) = nt + + @test max_remapped <= maximum(field) + @test min_remapped >= minimum(field) + + # Delete testing directory and files + rm(REGRID_DIR; recursive = true, force = true) + end + + @testset "test land_sea_mask for FT=$FT" begin + # Test setup + R = FT(6371e3) + test_space = TestHelper.create_space(FT, R = R) + comms_ctx = ClimaComms.SingletonCommsContext() + ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing + mkpath(REGRID_DIR) + + # Initialize dataset of all ones + data_path = joinpath(REGRID_DIR, "ls_mask_data.nc") + varname = "test_data" + TestHelper.gen_ncdata(FT, data_path, varname, FT(1)) + + # Test monotone masking + land_mask_mono = + Regridder.land_sea_mask(FT, REGRID_DIR, comms_ctx, data_path, varname, test_space, mono = true) + + # Test no new extrema are introduced in monotone remapping + nt = NCDataset(data_path) do ds + max_val = maximum(ds[varname]) + min_val = minimum(ds[varname]) + (; max_val, min_val) + end + (; max_val, min_val) = nt + + @test maximum(land_mask_mono) <= max_val + @test minimum(land_mask_mono) >= min_val + + # Test that monotone remapping a dataset of all ones conserves surface area + @test sum(land_mask_mono) - 4 * π * (R^2) < 10e-14 + + # Delete testing directory and files + rm(REGRID_DIR; recursive = true, force = true) + + # Set up testing directory + ispath(REGRID_DIR) ? rm(REGRID_DIR; recursive = true, force = true) : nothing + mkpath(REGRID_DIR) + + # Initialize dataset of all 0.5s + data_path = joinpath(REGRID_DIR, "ls_mask_data.nc") + varname = "test_data_halves" + TestHelper.gen_ncdata(FT, data_path, varname, FT(0.5)) + + # Test non-monotone masking + land_mask_halves = + Regridder.land_sea_mask(FT, REGRID_DIR, comms_ctx, data_path, varname, test_space, mono = false) + + # Masking of values below threshold should result in 0 + @test all(parent(land_mask_halves) .== FT(0)) + + # Delete testing directory and files + rm(REGRID_DIR; recursive = true, force = true) + end end end diff --git a/test/runtests.jl b/test/runtests.jl index cbd8bee10f..23fbb8c746 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,3 +3,4 @@ include("CouplerState/cplstate_interface.jl") # include("CoupledSimulations/cplsolver.jl") include("regridder_tests.jl") include("conservation_checker_tests.jl") +include("bcreader_tests.jl")