works without diffusion?

CliMA · Nov 21, 2023 · 59a39cf · 59a39cf
1 parent 02112d1
commit 59a39cf
Show file tree

Hide file tree

Showing 4 changed files with 65 additions and 53 deletions.
diff --git a/examples/freely_decaying_mediterraneum.jl b/examples/freely_decaying_mediterraneum.jl
@@ -1,46 +1,32 @@
 using GLMakie
 using Oceananigans
 using Oceananigans: architecture
+using Oceananigans.Fields: interpolate!
 using ClimaOcean
 using ClimaOcean.ECCO2
 using ClimaOcean.InitialConditions: three_dimensional_regrid!, adjust_tracers!
+using Oceananigans.ImmersedBoundaries: mask_immersed_field!
 using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
 using Oceananigans.Coriolis: ActiveCellEnstrophyConserving
 using Oceananigans.Units
+using Printf
 
 #####
 ##### Regional Mediterranean grid 
 #####
 
-function regrid_ecco_tracers(grid)
-    Tecco = ecco2_field(:temperature)
-    Secco = ecco2_field(:salinity)
-
-    ecco_tracers = (; Tecco, Secco)
-
-    # Make sure all values are extended properly before regridding
-    adjust_tracers!(ecco_tracers; mask = ecco2_center_mask(architecture(grid)))
-
-    T = CenterField(grid)
-    S = CenterField(grid)
-
-    # Regrid to our grid!
-    three_dimensional_regrid!(T, Tecco)
-    three_dimensional_regrid!(S, Secco)
-
-    return T, S
-end
-
 # A stretched vertical grid with a Δz of 1.5 meters in the first 50 meters
-z = stretched_vertical_faces(minimum_depth = 5000, 
-                             surface_layer_Δz = 1.75, 
+z = stretched_vertical_faces(depth = 5000, 
+                             surface_layer_Δz = 2.5, 
                              stretching = PowerLawStretching(1.070), 
                              surface_layer_height = 50)
 
-Nx = 20 * 42 # 1 / 20th of a degree
-Ny = 20 * 15 # 1 / 20th of a degree
+Nx = 4 * 42 # 1 / 4th of a degree
+Ny = 4 * 15 # 1 / 4th of a degree
 Nz = length(z) - 1
 
+@info "grid size: ($Nx, $Ny, $Nz)"
+
 grid = LatitudeLongitudeGrid(CPU();
                              size = (Nx, Ny, Nz),
                              latitude = (30, 45),
@@ -52,16 +38,23 @@ h = regrid_bathymetry(grid, height_above_water=1)
 
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(h))
 
-T, S = regrid_ecco_tracers(grid)
+Tecco, Secco = initial_ecco_tracers(architecture(grid))
+
+T = CenterField(grid)
+S = CenterField(grid)
+
+# Regrid to our grid!
+three_dimensional_regrid!(T, Tecco)
+three_dimensional_regrid!(S, Secco)
 
 mask_immersed_field!(T)
 mask_immersed_field!(S)
 
-# fig = Figure()
-# ax  = Axis(fig[1, 1])
-# heatmap!(ax, interior(T, :, :, Nz), colorrange = (10, 20), colormap = :thermal)
-# ax  = Axis(fig[1, 2])
-# heatmap!(ax, interior(S, :, :, Nz), colorrange = (35, 40), colormap = :haline)
+fig = Figure()
+ax  = Axis(fig[1, 1])
+heatmap!(ax, interior(T, :, :, Nz), colorrange = (10, 20), colormap = :thermal)
+ax  = Axis(fig[1, 2])
+heatmap!(ax, interior(S, :, :, Nz), colorrange = (35, 40), colormap = :haline)
 
 # Correct oceananigans
 import Oceananigans.Advection: nothing_to_default
@@ -78,11 +71,9 @@ model = HydrostaticFreeSurfaceModel(; grid,
                                       tracers  = (:T, :S, :e),
                                       coriolis = HydrostaticSphericalCoriolis(scheme = ActiveCellEnstrophyConserving()))
 
-set!(model, T = T, S = S)
-
-# Probably we'll need to diffuse? Probably not let's see now
+set!(model, T = T, S = S, e = 1e-6)
 
-simulation = Simulation(model, Δt = 20, stop_time = 2days)
+simulation = Simulation(model, Δt = 20, stop_iteration = 100)
 
 function progress(sim) 
     u, v, w = sim.model.velocities  
@@ -96,16 +87,16 @@ function progress(sim)
                    maximum(abs, T), maximum(abs, S))
 end
 
-simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))
+simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
 
-# warm up!
+# warm up for 100 iterations!
 run!(simulation)
 
 simulation.stop_time = 10*365days
 
-wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 2minutes, max_change = 1.1)
+wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1)
 
-simulation.callbacks = Callback(wizard, IterationInterval(10))
+simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
 
 simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                                                               indices = (:, :, Nz),

diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl
@@ -86,15 +86,9 @@ function regrid_bathymetry(target_grid;
     close(dataset)
 
     # Diagnose target grid information
-    Nxt, Nyt, Nzt = size(target_grid)
     arch = architecture(target_grid)
-    λt, φt, zt = nodes(target_grid, Face(), Face(), Face())
-
-    λ₁ = λt[1]
-    λ₂ = λt[Nxt]
-
-    φ₁ = φt[1]
-    φ₂ = φt[Nyt]
+    φ₁, φ₂ = y_domain(target_grid)
+    λ₁, λ₂ = x_domain(target_grid)
 
     # Calculate limiting indices on the bathymetry grid
     i₁ = searchsortedfirst(λ_data, λ₁)
@@ -151,24 +145,22 @@ function regrid_bathymetry(target_grid;
                                         latitude = (φ₁, φ₂),
                                         longitude = (λ₁, λ₂),
                                         z = (0, 1),
-                                        halo = halo_size(target_grid))
+                                        halo = (10, 10, 1))
 
     native_h = Field{Center, Center, Nothing}(native_grid)
     set!(native_h, h_data)
 
-    # 
     target_h = interpolate_bathymetry_in_passes(native_h, target_grid; passes = interpolation_passes)
 
     return target_h
 end
 
 # Here we can either use `regrid!` (three dimensional version) or `interpolate`
 function interpolate_bathymetry_in_passes(native_h, target_grid; passes = 10)
-
     Nλt, Nφt = Nt = size(target_grid)
     Nλn, Nφn = Nn = size(native_h)
 
-    if any(Nt .> Nn) # We are refining the grid (at least in one direction), more passes will not help!
+    if any(Nt[1:2] .> Nn[1:2]) # We are refining the grid (at least in one direction), more passes will not help!
         new_h = Field{Center, Center, Nothing}(target_grid)
         interpolate!(new_h, native_h)
         return new_h
@@ -177,6 +169,9 @@ function interpolate_bathymetry_in_passes(native_h, target_grid; passes = 10)
     latitude  = y_domain(target_grid)
     longitude = x_domain(target_grid)
 
+    @show latitude, longitude
+    @show x_domain(native_h.grid), y_domain(native_h.grid)
+
     ΔNλ = floor((Nλn - Nλt) / passes)
     ΔNφ = floor((Nφn - Nφt) / passes)
 
@@ -191,6 +186,7 @@ function interpolate_bathymetry_in_passes(native_h, target_grid; passes = 10)
 
     for pass = 1:passes - 1
         new_size = (Nλ[pass], Nφ[pass], 1)
+        @show Nt, Nn, new_size
 
         @info "pass number $pass with size $new_size"
         new_grid = LatitudeLongitudeGrid(size = new_size, 

diff --git a/src/DataWrangling/ECCO2.jl b/src/DataWrangling/ECCO2.jl
@@ -1,6 +1,8 @@
 module ECCO2
 
-export ecco2_field, ecco2_center_mask
+export ecco2_field, ecco2_center_mask, initial_ecco_tracers
+
+using ClimaOcean.InitialConditions: adjust_tracers!
 
 using Oceananigans
 using Oceananigans.BoundaryConditions
@@ -152,5 +154,28 @@ function ecco2_bottom_height_from_temperature()
     return bottom_height
 end
 
+function initial_ecco_tracers(architecture; 
+                              overwrite_existing = true, 
+                              initial_condition_file = "../data/initial_ecco_tracers.nc")
+
+    if overwrite_existing || !isfile(initial_condition_file)
+        T = ecco2_field(:temperature; architecture)
+        S = ecco2_field(:salinity; architecture)
+
+        # Make sure all values are extended properly before regridding
+        adjust_tracers!((; T, S); mask = ecco2_center_mask(architecture))
+
+        nc = Dataset(initial_condition_file, "w")
+        nc["T"] = interior(Tecco)
+        nc["S"] = interior(Tecco)
+    else 
+        nc = Dataset(initial_condition_file)
+        T = nc["T"]
+        S = nc["S"]
+    end
+
+    return T, S
+end
+
 end # module
 
diff --git a/src/InitialConditions.jl b/src/InitialConditions.jl
@@ -227,19 +227,19 @@ function three_dimensional_regrid!(a, b)
     source_size = Ns = size(source_grid)
 
     # Start by regridding in z
-    @info "Regridding in z"
+    @debug "Regridding in z"
     zgrid   = construct_grid(typeof(target_grid), arch, (Ns[1], Ns[2], Nt[3]), (xs, ys, zt), topo)
     field_z = Field(location(b), zgrid)
     regrid!(field_z, zgrid, source_grid, b)
 
     # regrid in y 
-    @info "Regridding in y"
+    @debug "Regridding in y"
     ygrid   = construct_grid(typeof(target_grid), arch, (Ns[1], Nt[2], Nt[3]), (xs, yt, zt), topo)
     field_y = Field(location(b), ygrid);
     regrid!(field_y, ygrid, zgrid, field_z);
 
     # Finally regrid in x
-    @info "Regridding in x"
+    @debug "Regridding in x"
     regrid!(a, target_grid, ygrid, field_y)
 end