Merge branch 'ss/new-zstar' of github.com:CliMA/Oceananigans.jl into …

…ss/new-zstar

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.95.0"
+version = "0.95.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/BoundaryConditions/polar_boundary_condition.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Grids: inactive_node, new_data
+using Oceananigans.Grids: inactive_node, new_data, YFlatGrid
 using CUDA: @allowscalar
 
 struct PolarValue{D, S}
@@ -20,6 +20,10 @@ end
 # Just a column
 @inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[1, 1, k]
 
+# YFlat grids do not have boundary conditions!
+latitude_north_auxiliary_bc(::YFlatGrid, args...) = nothing
+latitude_south_auxiliary_bc(::YFlatGrid, args...) = nothing
+
 # TODO: vectors should have a different treatment since vector components should account for the frame of reference.
 # For the moment, the `PolarBoundaryConditions` is implemented only for fields that have `loc[1] == loc[2] == Center()`, which
 # we assume are not components of horizontal vectors that would require rotation. (The `w` velocity if not a tracer, but it does

src/BuoyancyFormulations/BuoyancyFormulations.jl

@@ -1,12 +1,9 @@
 module BuoyancyFormulations
 
 export
-    BuoyancyForce, BuoyancyTracer, SeawaterBuoyancy, buoyancy_perturbationᶜᶜᶜ,
-    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState, TEOS10,
+    BuoyancyForce, BuoyancyTracer, SeawaterBuoyancy, LinearEquationOfState,
     ∂x_b, ∂y_b, ∂z_b, buoyancy_perturbationᶜᶜᶜ, x_dot_g_bᶠᶜᶜ, y_dot_g_bᶜᶠᶜ, z_dot_g_bᶜᶜᶠ,
-    top_buoyancy_flux,
-    buoyancy_frequency_squared,
-    BuoyancyField
+    top_buoyancy_flux, buoyancy, buoyancy_frequency, BuoyancyField
 
 using Printf
 using Oceananigans.Grids

src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

@@ -63,6 +63,7 @@ default_free_surface(grid; gravitational_acceleration=g_Earth) =
                                 free_surface = default_free_surface(grid, gravitational_acceleration=g_Earth),
                                 forcing::NamedTuple = NamedTuple(),
                                 closure = nothing,
+                                timestepper = :QuasiAdamsBashforth2,
                                 boundary_conditions::NamedTuple = NamedTuple(),
                                 tracers = (:T, :S),
                                 particles::ParticlesOrNothing = nothing,
@@ -93,6 +94,8 @@ Keyword arguments
                preallocated `CenterField`s.
   - `forcing`: `NamedTuple` of user-defined forcing functions that contribute to solution tendencies.
   - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
+  - `timestepper`: A symbol that specifies the time-stepping method. 
+                   Either `:QuasiAdamsBashforth2` (default) or `:SplitRungeKutta3`.
   - `boundary_conditions`: `NamedTuple` containing field boundary conditions.
   - `particles`: Lagrangian particles to be advected with the flow. Default: `nothing`.
   - `biogeochemistry`: Biogeochemical model for `tracers`.

src/OutputReaders/field_time_series.jl

@@ -13,7 +13,8 @@ using Oceananigans.Architectures
 using Oceananigans.Grids
 using Oceananigans.Fields
 
-using Oceananigans.Grids: topology, total_size, interior_parent_indices, parent_index_range
+using Oceananigans.Grids: topology, total_size, interior_parent_indices, parent_index_range, AbstractGrid
+using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
 
 using Oceananigans.Fields: interior_view_indices, index_binary_search,
                            indices_summary, boundary_conditions
@@ -386,6 +387,10 @@ function FieldTimeSeries(loc, grid, times=();
                                        times, path, name, time_indexing, reader_kw)
 end
 
+isreconstructed(grid::JLD2.ReconstructedStatic) = true
+isreconstructed(grid::AbstractGrid) = false
+isreconstructed(grid::ImmersedBoundaryGrid) = isreconstructed(grid.underlying_grid)
+
 """
     FieldTimeSeries{LX, LY, LZ}(grid::AbstractGrid [, times=()]; kwargs...)
 
@@ -466,8 +471,6 @@ function FieldTimeSeries(path::String, name::String;
         (:, :, :)
     end
 
-    isnothing(grid) && (grid = file["serialized/grid"])
-
     if isnothing(architecture) # determine architecture
         if isnothing(grid) # go to default
             architecture = CPU()
@@ -476,11 +479,54 @@ function FieldTimeSeries(path::String, name::String;
         end
     end
 
-    if boundary_conditions isa UnspecifiedBoundaryConditions
-        boundary_conditions = file["timeseries/$name/serialized/boundary_conditions"]
-        boundary_conditions = on_architecture(architecture, boundary_conditions)
-    end
+    isnothing(grid) && (grid = file["serialized/grid"])
 
+    # If isreconstructed(grid), it probably means that the data was generated prior to
+    # Oceananigans version 0.95.0 (12/13/2024). In this case, the best we can do is to try to rebuild
+    # the grids manually using the non-serialized grid data. Here, we support RectilinearGrid
+    # and LatitudeLongitudeGrid (but not OrthogonalSphericalShellGrid) and we also assume
+    # GridFittedBottom if the grid is an ImmersedBoundaryGrid. If these assumptions can be relaxed
+    # in the future, they should.
+    if isreconstructed(grid)
+        isibg = grid isa ImmersedBoundaryGrid
+        test_grid = isibg ? grid.underlying_grid : grid
+        address = isibg ? "grid/underlying_grid" : "grid"
+        Nx = file["$address/Nx"]
+        Ny = file["$address/Ny"]
+        Nz = file["$address/Nz"]
+        Hx = file["$address/Hx"]
+        Hy = file["$address/Hy"]
+        Hz = file["$address/Hz"]
+        zᵃᵃᶠ = file["$address/zᵃᵃᶠ"]
+        z = file["$address/Δzᵃᵃᶠ"] isa Number ? (zᵃᵃᶠ[1], zᵃᵃᶠ[Nz+1]) : zᵃᵃᶠ[1:Nz+1]
+        topo = topology(grid)
+        size = (Nx, Ny, Nz)
+        halo = (Hx, Hy, Hz)
+
+        if :λᶜᵃᵃ ∈ propertynames(test_grid) # I guess its a LatitudeLongitudeGrid.
+            λᶠᵃᵃ = file["$address/λᶠᵃᵃ"]
+            φᵃᶠᵃ = file["$address/φᵃᶠᵃ"]
+            λ = file["$address/Δλᶠᵃᵃ"] isa Number ? (λᶠᵃᵃ[1], λᶠᵃᵃ[Nx+1]) : λᶠᵃᵃ[1:Nx+1]
+            φ = file["$address/Δφᵃᶠᵃ"] isa Number ? (φᵃᶠᵃ[1], φᵃᶠᵃ[Ny+1]) : φᵃᶠᵃ[1:Ny+1]
+            domain = (latitude=φ, longitude=λ, z=z)
+            underlying_grid = LatitudeLongitudeGrid(architecture; size, halo, topology=topo, domain...)
+        else
+            xᶠᵃᵃ = file["$address/xᶠᵃᵃ"]
+            yᵃᶠᵃ = file["$address/yᵃᶠᵃ"]
+            x = file["$address/Δxᶠᵃᵃ"] isa Number ? (xᶠᵃᵃ[1], xᶠᵃᵃ[Nx+1]) : xᶠᵃᵃ[1:Nx+1]
+            y = file["$address/Δyᵃᶠᵃ"] isa Number ? (yᵃᶠᵃ[1], yᵃᶠᵃ[Ny+1]) : yᵃᶠᵃ[1:Ny+1]
+            domain = (; x, y, z)
+            underlying_grid = RectilinearGrid(architecture; size, halo, topology=topo, domain...)
+        end
+
+        if isibg
+            bottom_height = file["grid/immersed_boundary/bottom_height"]
+            bottom_height = view(bottom_height, 1+Hx:Nx+Hx, 1+Hy:Ny+Hy, 1)
+            grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height))
+        else
+            grid = underlying_grid
+        end
+    end
 
     # This should be removed eventually... (4/5/2022)
     grid = try
@@ -529,6 +575,11 @@ function FieldTimeSeries(path::String, name::String;
         end
     end
 
+    if boundary_conditions isa UnspecifiedBoundaryConditions
+        boundary_conditions = file["timeseries/$name/serialized/boundary_conditions"]
+        boundary_conditions = on_architecture(architecture, boundary_conditions)
+    end
+
     close(file)
 
     LX, LY, LZ = Location

src/OutputReaders/set_field_time_series.jl

@@ -36,6 +36,7 @@ function set!(fts::InMemoryFTS, path::String=fts.path, name::String=fts.name)
 
         # Note: use the CPU for this step
         field_n = Field(location(fts), path, name, file_iter,
+                        grid = on_architecture(CPU(), fts.grid),    
                         architecture = cpu_architecture(arch),
                         indices = fts.indices,
                         boundary_conditions = fts.boundary_conditions)