Always interpolate bathymetry in passes (#496)

* Always interpolate bathymetry in passes

* Update test_bathymetry.jl

* Update test_ocean_sea_ice_model.jl

---------

Co-authored-by: Simone Silvestri <[email protected]>

Author: glwagner

src/Bathymetry.jl

@@ -224,38 +224,6 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
     Nλt, Nφt = Nt = size(target_grid)
     Nλn, Nφn = Nn = size(native_z)
 
-    resxt = minimum_xspacing(target_grid)
-    resyt = minimum_yspacing(target_grid)
-
-    resxn = minimum_xspacing(native_z.grid)
-    resyn = minimum_yspacing(native_z.grid)
-
-    # Check whether we are refining the grid in any directions.
-    # If so, skip interpolation passes, as they are not needed.
-    if resxt < resxn || resyt < resyn
-        target_z = Field{Center, Center, Nothing}(target_grid)
-
-        @info "Interpolating bathymetry of size $Nn onto a $gridtype target grid of size $Nt"
-        interpolate!(target_z, native_z)
-
-        if passes > 0
-            passes_str = passes == 1 ? "1 interpolation pass" : "$passes interpolation passes"
-
-            @info string("Skipping extra $passes_str for bathymetry.", '\n',
-                        "Extra interpolation passes may only be used to coarsen bathymetry", '\n',
-                        "and require that the bathymetry is finer than the target grid in", '\n',
-                        "both horizontal directions.", '\n',
-                        "Grid info:", '\n',
-                        "target grid", '\n',
-                        "├── minimum x-spacings ", @sprintf("%.3e", resxt), " m", '\n',
-                        "└── minimum y-spacings ", @sprintf("%.3e", resyt), " m", '\n',
-                        "bathymetry grid", '\n',
-                        "├── minimum x-spacings ", @sprintf("%.3e", resxn), " m", '\n',
-                        "└── minimum y-spacings ", @sprintf("%.3e", resyn), " m")
-        end
-        return target_z
-    end
-
     # Interpolate in passes
     latitude  = y_domain(native_z.grid)
     longitude = x_domain(native_z.grid)
@@ -275,7 +243,7 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
     @info "Interpolation passes of bathymetry size $(size(old_z)) onto a $gridtype target grid of size $Nt:"
     for pass = 1:passes - 1
         new_size = (Nλ[pass], Nφ[pass], 1)
-        @info "    pass $pass to size $new_size"
+        @info "    Performing pass $pass to size $new_size."
 
         new_grid = LatitudeLongitudeGrid(architecture(target_grid), Float32,
                                          size = new_size,
@@ -469,4 +437,4 @@ end
 retrieve_bathymetry(grid, ::Nothing; kw...) = regrid_bathymetry(grid; kw...)
 retrieve_bathymetry(grid; kw...)            = regrid_bathymetry(grid; kw...)
 
-end # module
+end # module

test/test_bathymetry.jl

@@ -6,8 +6,8 @@ using ClimaOcean
 
 using ClimaOcean.Bathymetry: remove_minor_basins!
 
-@testset "Availability of Bathymetry" begin
-    @info "Testing Bathymetry utils..."
+@testset "Bathymetry construction and smoothing" begin
+    @info "Testing Bathymetry construction and smoothing..."
     for arch in test_architectures
 
         grid = LatitudeLongitudeGrid(arch;
@@ -50,18 +50,6 @@ using ClimaOcean.Bathymetry: remove_minor_basins!
         # of -1000 meters
         @test mean(view(bottom_height, 51:100, :, 1)) == -1000
 
-        grid = LatitudeLongitudeGrid(arch;
-                                     size = (200, 200, 10),
-                                     longitude = (0, 2),
-                                     latitude = (-0.1, 0.1),
-                                     z = (-6000, 0))
-
-        control_bottom_height = regrid_bathymetry(grid)
-        interpolated_bottom_height = regrid_bathymetry(grid; interpolation_passes=100)
-
-        # Testing that multiple passes do not change the solution when refining the grid
-        @test parent(control_bottom_height) == parent(interpolated_bottom_height)
-
         grid = LatitudeLongitudeGrid(arch;
                                      size = (200, 200, 10),
                                      longitude = (0, 100),

test/test_ocean_sea_ice_model.jl

@@ -58,7 +58,7 @@ using ClimaSeaIce.Rheologies
 
             bottom_height = regrid_bathymetry(grid;
                                               minimum_depth = 10,
-                                              interpolation_passes = 20,
+                                              interpolation_passes = 5,
                                               major_basins = 1)
 
             grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)